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SUMMARY 


In  a  general  sense,  this  report  represents  an  effort  to  clarify  the  relationship  of 
discrete  and  continuous  wavelet  transforms.  More  narrowly,  it  focuses  on  bringing 
together  two  separately  motivated  implementations  of  the  wavelet  transform,  the 
algorithme  'a  trous  and  Mallat’s  multiresolution  decomposition.  It  is  observed  that 
these  algorithms  are  both  special  cases  of  a  single  filter  bank  structure,  the  discrete 
wavelet  transform,  the  behavior  of  which  is  governed  by  one’s  choice  of  filters.  In 
fact,  the  'a  trous  algorithm,  originally  devised  as  a  computationally  efficient  imple¬ 
mentation,  is  more  properly  viewed  as  a  nonorthogonal  multiresolution  algorithm  for 
which  the  discrete  wavelet  transform  is  exact.  Moreover,  it  is  shown  that  the  com¬ 
monly  used  Lagrange  'a  trous  filters  are  in  one-to-one  correspondence  with  the  convo¬ 
lutional  squares  of  the  Daubechies  filters  for  orthonormal  wavelets  of  compact  sup¬ 
port. 

A  systematic  framework  for  the  discrete  wavelet  transform  is  provided,  and  con¬ 
ditions  are  derived  under  which  it  computes  the  continuous  wavelet  transform  exactly. 
Suitable  filter  contraints  for  finite  energy  and  boundedness  of  the  discrete  transform 
are  also  derived.  Finally,  relevant  signal-processing  parameters  are  examined,  and  it 
is  remarked  that  orthonormality  is  balanced  by  restrictions  on  resolution. 
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1.  INTRODUCTION 


Wavelets  are  rapidly  finding  application  as  a  tool  for  the  analysis  of  nonstation¬ 
ary  signals.  [1]  —  [5]  However,  with  the  notable  exception  of  orthonormal  wavelets 
[6]  —  [9],  very  little  literature  has  been  devoted  to  linking  discrete  implementations  to 
the  continuous  transform.  As  in  the  case  of  the  discrete  Fourier  transform,  these 
implementations  (or  filter  banks)  can,  and  should,  stand  on  their  own  as  abstract 
decompositions  of  discrete  time  series.  Their  wide  sweeping  significance,  however, 
lies  in  their  interpretation  as  wavelet  transforms.  In  a  general  sense,  this  report 
represents  an  effort  to  clarify  the  relationship  of  discrete  and  continuous  wavelet 
transforms.  More  narrowly,  it  focuses  on  bringing  together  two  separately  motivated 
implementations  of  the  wavelet  transform.  One  of  them,  the  algorithme  'a  trous 1  for 
nonorthogonal  wavelets,  [4]  —  [5],  was  developed  for  music  synthesis  [2]  and  is  partic¬ 
ularly  suitable  for  signal  processing.  The  other,  the  multiresolution  approach  of  S. 
Mallat  and  Y.  Meyer,  originally  used  in  Image  processing,  employs  orthonormal 
wavelets.  [6]  —  [10]  The  latter  algorithm,  apart  from  its  wavelet  interpretation,  was 
discovered  previously  in  the  form  of  quadrature  mirror  filter  (QMF)  filter  banks  with 
perfect  reconstruction  where  it  finds  application  in  speech  transmission  and  split-band 
coding.  [11]  —  [13] 

A  glance  at  these  two  algorithms  suffices  to  reveal  closely  related  structures.  In 
fact,  apart  from  the  constraints  on  their  filters,  the  decimated  'a  trous  [5]  and  Mallat 
algorithms  are  identical.  We  are  thus  led  to  examine  the  expanded  family  of  algo¬ 
rithms  encompassing  both  types  of  filters.  In  this  vein,  it  is  showp,  that  the  Lagrange 
interpolation  filters  commonly  employed  by  the  'a  trous  algorithm  are  actually  the 
squares  (in  a  convolutional  sense)  of  the  Daubechies  filters  for  compact  orthonormal 
wavelets.  We  also  derive  conditions  under  which  the  discrete  implementation  com¬ 
putes  a  continuous  wavelet  transform  exactly  and  find  that  they  bear  an  intimate  rela¬ 
tionship  to  the  'a  trous  constraints. 

From  a  more  general  viewpoint  the  situation  is  as  follows:  The  algorithms  to  be 
discussed  all  are  filter  bank  structures  (see  figure  1.1).  Their  only  distinguishing 
feature  is  the  choice  of  two  finite  length  filters,  a  lowpass  filter  f  and  a  bandpass  filter 
g.  The  lowpass  condition,  expressed  more  precisely  as  S  fk  =  V2,  is  necessary  to  the 
construction  of  a  corresponding  continuous  wavelet  function.  The  bandpass  require¬ 
ment,  while  apparently  not  essential  to  all  applications,  insures  that  finite  energy  sig¬ 
nals  lead  to  finite  energy  transforms  (cf.,  section  6).  Under  these  conditions  the  filter 
bank  output  will  be  referred  to  as  the  Discrete  Wavelet  Transform  (DWT),  a  termi¬ 
nology  which  will  become  clear  in  the  course  of  the  report. 

One  class  of  DWT  filter  pairs  are  the  Daubechies  filters  [8]  which  yield  orthogo¬ 
nal  wavelet  decompositions  and  constitute,  in  more  conventional  terms,  a  QMF  filter 
bank  with  perfect  reconstruction.  Another  is  that  for  which  the  lowpass  filter  satis¬ 
fies  the  ' a  trous  condition  f2k  =  <5(k)  /V2.  Such  filters,  which  simply  serve  to  interpo¬ 
late  every  other  point,  correspond  to  a  nonorthogonal  wavelet  decomposition.  As 
mentioned  above,  if  they  are  further  restricted  to  be  Lagrangian  interpolators,  they 

Literally,  "algorithm  with  holes”,  this  terminology,  taken  from  [5],  refers  to  the  fact  that  all 
the  even  coefficients  of  the  relevant  filter  (with  the  exception  of  the  center)  are  zero. 
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Figure  1.1  A  wavelet  filter  bank  structure.  The  down-arrow 
indicates  decimation.  The  output  of  the  transform  is  the 
family  of  signals  w1,  forming  the  two  parameter  transform  w„ 
in  the  scale-time  plane.  Following  terminology  to  be  intro¬ 
duced,  w1  is  the  (decimated)  discrete  wavelet  transform. 


become  the  squares  of  the  Daubechies  filters,  which  is  quite  remarkable  in  considera¬ 
tion  of  the  totally  different  derivations.  This  also  implies  that  a  maximally  flat  filter 
with  the  same  number  of  vanishing  derivatives  at  0  and  7r  is  a  Lagrangian  interpolator; 
a  fact  which  may  aid  in  the  design  of  maximally  flat  filters.  [14] 

A  fundamental  question  is  when  do  these  discrete  implementations  yield  exact 
(i.e.,  sampled)  versions  of  a  continuous  wavelet  transform?  Aside  from  regularity 
conditions  relating  to  smoothness  [8],  we  find  that  if  f  is  'a  trous,  then  the  DWT  coin¬ 
cides  with  a  continuous  wavelet  transform  by  a  wavelet  g(t)  whose  samples  g(n)  form 
the  filter  g.  Even  if  f  is  not  'a  trous,  the  algorithm  is  exact  provided  the  signal  lies  in 
an  appropriate  subspace;  however,  in  that  instance,  the  corresponding  wavelet 
depends  on  f  as  well  as  g.  This  is  the  situation  in  the  orthonormal  case  where,  more¬ 
over,  the  filter  g  is  almost  completely  determined  from  f  through  the  constraints  of 
orthogonality. 

In  the  remainder  of  this  introduction  we  define,  and  briefly  motivate,  wavelet 
transforms  at  various  levels  of  discretization.  Section  2  contains  an  abbreviated 
derivation  of  the  'a  trous  algorithm  followed  by  a  description  of  the  Mallat  algorithm. 
(The  uninitiated  reader  is  particularly  referred  to  [1],  [6],  [8].)  In  section  3  we  define 
the  undecimated  DWT,  relate  it  to  the  decimated  transform,  and  provide  algorithms 
for  its  computation.  Section  4  states  and  proves  several  theorems  which  delineate  the 
relationship  between  the  DWT  and  the  continuous  wavelet  transform.  It  may  be  read 
independently  of  the  algorithms  of  section  2  although  the  motivation  for  the  construc¬ 
tions  may  not  be  clear.  Section  5  defines  the  Lagrange  'a  trous  filters  and  proves  that 
they  are  the  squares  of  the  Daubechies  filters.  In  section  6,  we  formulate  the  inver¬ 
sion  problem  and  provide  filter  constraints  which  insure  finite  energy  and  bounded 
operators.  It  concludes  with  a  short  examination  of  the  tradeoffs  involved  in  choosing 
the  bandpass  filter,  emphasizing  the  differences  of  the  orthonormal  and  nonorthogo- 
nal  cases. 

TRANSFORM  DEFINITIONS 

The  continuous  wavelet  transform  of  a  signal  s(t)  takes  the  form 
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(1.1) 


=  VT  /  8<i7A)  s(t)  dl 

where  g  is  the  analyzing  wavelet,  a  represents  a  time  dilation,  b  a  time  translation, 
and  the  bar  stands  for  complex  conjugate.  The  normalization  factor  1/Va  is  perhaps 
most  effectively  visualized  as  endowing  |  W (a,b)  |2  with  units  of  power/Hz. z  Certain 
weak  "admissibility"  conditions  are  usually  required  on  g(t)  for  it  to  be  a  candidate  for 
an  analyzing  wavelet;  namely,  square  integrability  and 

/  da;  <  oo  (1.2) 

M 

where  g(w)  is  the  Fourier  transform  of  g(t).  They  insure  that  the  transformation  is  a 
bounded  invertible  operator  in  the  appropriate  spaces  (cf.,  J8],  [16]).  If  g(cu)  is  dif¬ 
ferentiable,  then  it  suffices  that  g  be  zero  mean,  i.e.,  that  J  g(t)dt  =0,  for  equation 
(1.2)  to  be  satisfied. 

In  the  realm  of  signal  processing,  the  significance  of  equation  (1.1)  is  probably 
(or,  at  least,  traditionally)  best  grasped  by  comparing  it  to  the  short-time  Fourier 
transform: 


F(u,  b)  =  /  h(t  -  b)  eiw‘  s(t)  dt  .  (1.3) 

Thus,  to  obtain  F(w,  b),  one  multiplies  the  signal  by  an  appropriate  window  h  (such 
as  a  Gaussian)  centered  at  time  b  and  then  takes  the  Fourier  transform.  In 
mathematical  terms,  equation  (1.3)  is  an  expansion  of  the  signal  in  terms  of  a  family 
of  functions  h(t  —  b)elut,  which  are  generated  from  a  single  function  h(t)  through 
translations  b  in  time  and  translations  u  in  frequency.  In  contrast,  the  wavelet 

transform  (1.1)  is  an  expansion  in  functions  g(- - — )  generated  by  translations  b  in 

...  3  ^ 

time  and  dilations  a  in  time.  Thus,  the  continuous  wavelet  transform  resembles  a 
(continuous)  bank  of  short-time  Fourier  transforms  with  a  different  window  for  each 
frequency.  The  significance  of  this  is  that,  while  the  basis  functions  in  (1.3)  have  the 
same  time  and  frequency  resolution  (that  of  h(t)  and  h(w))  at  all  points  of  the 
transform  plane,  those  of  (1.1)  have  time  resolution  (that  of  g(t/fl))  which  decreases 
with  a  and  frequency  resolution  (that  of  g (aw))  which  increases  with  a.  This  property 
can  be  a  great  advantage  in  signal  processing  since  high-frequency  signal 


“  The  energy  density  is  given  by  |  W  | :  y— .  an  expression  which  is  intimately  linked  to 

Q 

representations  of  the  affine  group  (cf.,  [3],  [15]  -  [16]).  Since  a  is  proportional  to 
bandwidth,  |  W  j2  is  power/Hz.  See  section  6  for  a  discussion  of  the  discrete  case. 

Alternatively,  dilations  in  time  may  be  considered  contractions  in  frequency  since  the 
Fourier  transform  of  g(t/a)  is  ag(<i  jj) .) 
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characteristics  are  generally  highly  localized  in  time  whereas  slowly  varying  signals 
require  good  low-frequency  resolution. 

As  originally  proposed  by  Morlet  et  al.  [17],  g  was  a  modulated  Gaussian 

g(t)  =  e1"0*  e-‘2/2  ,  (1.4) 

and  this  function  is  still  the  prototypical  analyzing  wavelet  for  signal-processing  appli- 

— (u>— — )2  a2/ 2 

cations.  [1]  The  window  function  g(t /a)  has  Fourier  transform  g(aw)  =ae  a  , 
which  has  analysis  frequency  t/0 /a.  We  emphasize  that  i/0  is  simply  a  parameter 
which  determines  the  analyzing  wavelet;  its  role  should  not  be  confused  with  that  of  a 
even  though  the  scale  axis  is  often  expressed  in  terms  of  frequency  under  the  transfor¬ 
mation  a  —*■  fa.  Observe  that,  (1.4)  only  satisfies  the  admissibility  condition  (1.2) 
approximately  (cf.,  [16],  [18]).  As  expected  its  bandwidth  is  proportional  to  1/a,  thus 
giving  rise  to  a  constant  relative  bandwidth;  i.e.,  BW /(v0/a)  =  constant.  This  feature 
is  also  reflected  in  the  narrowing  of  the  time  window  at  higher  frequencies;  i.e.,  at 
smaller  a.  In  general,  one’s  choice  of  the  function  g  is  dictated  by  its  time  and  fre¬ 
quency  localization  properties  (see  [3],  [18]). 

We  shall  be  exclusively  concerned  with  discrete  values  for  a  and  b.  In  particu¬ 
lar,  we  assume  that  a  has  the  form  a  =  21  where  i  is  termed  the  octave  of  the 
transform.  The  integral  (1.1)  then  yields  a  wavelet  series 

W(2S  n)  4  -±-  J  s(t)  dt  .  (1.5a) 


We  remark  that  finite  energy  for  the  wavelet  transform  is  not  at  all  equivalent  to  finite 
energy  for  the  wavelet  series.  It  depends  on  the  sampling  grid  as  well  as  the  function 
g(t).  [3]  Thus,  the  admissibility  condition  (1.2)  is  not  necessarily  appropriate  in  the 
discrete  case  and  shall  be  replaced  with  conditions  on  the  relevant  filters  in  section  6. 
In  addition,  we  shall  often  take  b  to  be  a  multiple  of  a4 

W(2\  2‘n)  4  -L-  /  g(i-  -  n)  s(t)  dt  .  (1.5b) 


A  logical  step  in  applying  the  theory  to  discrete  signals  is  to  discretize  the 
integial  in  (1.5) 


Physically,  this  reflects  a  need  for  less  frequent  sampling  of  the  transform  output  at  lower 
frequencies  (i.e.,  larger  scales  a).  Mathematically,  b  =-  2'n  has  its  roots  in  the  orthonormal 
wavelets  where  it  suffices  for  invertibility  of  the  transform.  [6]  The  general  case,  however,  is 
much  more  complex.  [3]  Too  sparse  a  sampling  leads  to  incompleteness;  oversampling  results 
in  a  redundant  set  of  functions. 
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w<2'.  2‘n)  s(k) 


(1.6) 


The  sample  rate  has  been  set  equal  to  one.  As  indicated  by  2’n  on  the  left  hand  side, 
(1.6),  as  well  as  (1.5b),  are  decimated  wavelet  transforms.  Octave  i  is  only  output 
every  2‘  samples.  In  this  form  the  resulting  algorithms  will  not  be  translation  invari¬ 
ant.  [7]  This  is  easily  seen  by  substituting  s(k—  r)  for  s(k)  which  produces 
w(2’,  2‘(n— r/21)),  an  integer  translation  of  w(2',  2'n)  only  if  r  is  a  multiple  of  2\ 
However,  the  invariance,  which  is  lost  by  decimation,  is  easily  restored  by  separately 
filtering  the  even  and  odd  sequences  (cf.,  section  3)  or  by  using  an  equivalent  algo¬ 
rithm,  also  described  in  section  3.  Our  major  reason  for  starting  with  (1.5b)  rather 
than  (1.5a)  is  historical.  It  delineates  the  relationship  of  the  DWT  to  the  QMF  filter 
banks  and  (orthonormal)  wavelet  structures  already  found  in  the  literature.  [6]  —  [9] 
It  also  simplifies  our  derivation  of  the  'a  trous  algorithm  and  readily  lends  itself  to 
physical  interpretation  (cf.,  section  6).  Note  that  the  original  Mallat  algorithm  [6]  was 
decimated;  ' a  trous  [4]  was  not. 

Proceeding  from  (1.6),  we  shall  be  able  to  to  arrive  at  the  DWT  pictured  in  fig¬ 
ure  1.1, 

[**]«.  ^  ?  f2n-j  [Si_1]j  (1-7) 

) 

M„  =  ?  gn-j  [Si_1]j  • 

J 

where  [w‘]n  corresponds  to  w(2J,  2Jn)  of  equation  (1.6)  and  s°  is  the  original  signal  s. 
The  mysterious  appearance  of  the  filter  f  in  (1.7)  will  be  unraveled  in  the  derivation  of 
the  a  trous  algorithm  in  section  2.  Finally,  we  shall  come  full  circle  in  section  4 
where,  under  quite  general  conditions,  we  show  the  existence  of  a  function  g(t)  with 
g(n)  =  g„  £  g_n  and  such  that  the  DWT  acting  on  the  sampled  signal  is  exactly  the 
sampled  output  of  the  continuous  wavelet  transform  (i.e.,  of  the  wavelet  series).5  In 
other  words,  the  DWT  with  filter  g  defined  by  g„  =  gA(n),  which  was  originally  con¬ 
ceived  as  an  approximation  of  the  (continuous)  WT  for  an  arbitrary  analyzing  wavelet 
gA(t),  is  exact  for  another  wavelet  function  gB(t)  where  gB(n)  =  g„  for  all  n.  Of 
course,  if  there  is  sufficient  regularity,  gA(t)  and  gB(t)  will  be  close  since  they  coincide 
on  the  integers  up  to  the  length  of  g. 

Before  embarking  on  this  voyage,  we  summarize,  and  hopefully  clarify,  the 
plethora  of  transforms  with  a  brief  analogy  to  the  Fourier  transform,  Fourier  series, 
the  discretized  z-transform,  and  the  discrete  Fourier  transform  (DFT).  The  Fourier 
transform  of  a  continuous  signal  s(t). 


SM 


oc 

£  f  s(t)  dt  , 

— 00 


(1.8) 


5  The  adjoint  filter  gk  -  g_k  is  used  to  simplify  future  notation.  It  corresponds  to  the  in¬ 
tegrand  g( — t)  *  s  found  in  equations  (1.1),  (1.5),  and  (1.6). 
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is  a  function  of  the  continuous  variable  u.  Restricting  it  to  a  discrete  (one¬ 
dimensional)  grid  results  in  the  coefficients  of  a  Fourier  series 


S(27rm)  =  /  e-2!,imt  s(t)  dt,  (1.9) 

—OO 

which  in  turn  may  be  computed  approximately  by 

sz(2;rmAt)  =  £  e2*imk*t  s(kAt)At,  (1.10) 

k 

the  z-transform  of  sn  s(nAt)  output  at  discrete  points  e2’rimAt.  If  s(t)  is  band-limited 
and  sampled  at  an  appropriate  rate,  At  =  1/N,  then  the  above  may  be  computed 
exactly  using  the  DFT 


N  2ffimk 

sm  =  35-  N  sk  .  (1.11) 

These  correspond  precisely  to  W (a,  b ),  W(2',  n),  w(2',  n),  and  undecimated  w'n.  With 
wavelets,  however,  we  have  the  additional  difficulty  of  dealing  with  a  whole  class  of 
functions  g(t)  rather  than  simply  e'wt.  Also  complicating  things  are  its  two- 
dimensional  structure  and  the  decimated  versions,  which,  due  to  their  2'n  dependency 
on  i,  play  a  distinguished  role  without  analogy  in  the  one-dimensional  case. 
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2.  TWO  ALGORITHMS 


NOTATION 

Decimation,  which  appears  as  a  down  arrow  in  figure  1.1,  plays  a  pivotal  role  in 
all  DWT  algorithms.  However,  it  leads  to  operators  which  are  not  time  invariant  and 
present  a  potential  source  of  confusion.  It  is  thus  worthwhile  to  first  establish  some 
formal  notation. 

Signals  and  filters  in  boldface  will  be  treated  as  vectors,  in  which  case  *  indi¬ 
cates  discrete^  convolution  and  yields  a  vector.  The  symbol  t  will  be  used  for  adjoint 
filter  (ff)k  =  f_k.  Note  that  this  is  the  complex  conjugate  reversal  and  does  not  imply 
a  conversion  of  a  column  vector  to  a  row  vector.  The  above  mentioned  decimation 
operator  is  represented  by  a  matrix 

Akj  4  6(2k  -  j) 

=  *2k>j  (2.1) 

where  <5;:  is  the  Kroneker  delta  and  <5(k)  £  S^q.  Also  of  significance  is  the  dilation 
operator 

Dkj  £  6(k  -  2j) 

=  Skt2j  ,  (2-2) 

which  dilates  a  vector  by  inserting  zeros.  Observe  that  A  and  D  are  transposes  of 
each  other,  and  that  although  they  are  linear,  they  are  not  time  invariant;  i.e.,  they 
are  not  functions  of  k  —  j. 

Convolution  followed  by  decimation  becomes  [A(f  *  s)]k  =  S  Akj  [f  *  s]j  = 
[f  *  skk  =  ^  f2k-msm-  However,  a  particularly  insidious  pitfall  remains,  namely, 


(Af  )  *  s  *  A  (f  *  s)  .  (2.3) 

This  associativity  problem  may  be  avoided  by  replacing  convolution  by  f  with  a 
matrix,  F,  defined  by 

F,,  A  f,_j  ,  (2.4) 

which  shall  occasionally  be  used  in  our  proofs.  A  trivial  calculation  yields 
AFs  =  A  (f  *  s).  The  symbol  t  will  also  be  usedjor  the^ adjoint  of  matrices.  This  is 
consistent  with  the  above  notation  where  (F+) ^  £  Fj;  =  fj_;  4  f;+_j. 

We  define  the  Fourier  transform  of  a  function  s(t)  by 
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s (u>)  4  /  s(t)  e  ,u,t  dt  , 


(2.5) 


and  the  z-transform  (on  the  unit  circle)  of  a  discrete  signal  s  by 


sz(-’)  =  Ssn 

n 


e-iun 


(2.6) 


In  the  subsequent  interplay  between  continuous  and  discrete  functions  one  must  be 
careful  to  distinguish  the  usage  of  these  two  transforms.  Ignoring  their  differences 
can  easily  lead  to  erroneous  conclusions.  In  particular, 

(Af)2M  =  E  f2„  e-‘"“  -  £  f„  *  f,(2w)  .  (2.7) 

n  n  even 


Some  comment  concerning  filter  definitions  is  also  appropriate.  Usage  in  the 
literature  is  uniform  only  up  to  the  adjoint.  Also,  the  z-transform  is  sometimes 
defined  with  a  positive  exponential  which  leads  to  similar  differences  in  the  frequency 
domain.  In  keeping  with  signal-processing  applications  we  have  chosen  (2.6)  as 
above,  consistent  with  the  Fourier  transform,  and  we  shall  define  our  filters  so  that 
adjoints  do  not  appear  in  convolutions.  This  produces  a  minimum  of  adjoints  and 
greatly  simplifies  the  notation.  Unfortunately,  it  also  results  in  the  definition 
Sn  =  g(n)  anc*  the  introduction  of  ft  as  an  interpolation  filter  whereas  g  and  f  would 
be  more  natural.  Note,  also,  that  our  filters  are  the  adjoints  of  the  filters  defined  by 
I.  Daubechies  [8],  although  their  z-transforms  coincide  since  [8]  defines  the  z- 
transform  with  a  plus  sign. 

THE  A  TROUS  ALGORITHM 

We  take  the  discretized  wavelet  series  (1.6)  as  our  starting  point.  The  difficulty 
in  implementing  (1.6)  is  that,  even  for  g(t)  of  finite  support,  as  i  increases,  g  must  be 
sampled  at  progressively  more  points,  creating  a  large  computational  burden.  The 
solution  posed  in  [4]  is  to  approximate  the  values  at  nonintegral  points  through  inter¬ 
polation  via  a  finite  filter  f'.  The  resulting  recursion  is  highly  efficient  and  may  be 
implemented  with  the  filter  bank  structure  of  figure  1.1. 

The  interpolation  is  perhaps  best  introduced  with  an  example.  Let  f+  be  the 
filter  (0.5,  1.0,  0.5).  Then, 


S  -  2k  g(k) 

k 


1  ,  /  n  1  \  .  /  n+1  x  \ 

y(  g(~  )  +  g(~  )  ) 


n  even 

n  odd 


(2.8) 


approximates  a  sampling  of  g(t/2).  With  the  help  of  the  dilation  operator  D,  this  may 
Again,  we  use  the  adjoint  in  our  definition  to  simplify  subsequent  notation. 
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be  formalized  as  a  general  procedure  for  dyadic  interpolation.  The  steps  are  illus¬ 
trated  in  figure  2.1.  Let  g  be  a  filter  defined  by  g*  4  g(n).  First,  we  spread  g*  to 
pro- 


g(n):  g(-l)  g(0)  g(l)  g(2) 


Dgt  : 

g(-i)  o 

g(0)  o 

gd) 

0  g(2) 

jfV 2 

»  g(n/2): 

g(-l)  x 

g(0)  x 

gd) 

x  g(2) 

Figure  2.1  Diagram  illustrating  the  dialation  and  interpolation  of  a 
function  g(t)  :  g(n/2)  ~  V2f*  *  (Dg^. 


vide  space  in  which  to  put  the  interpolated  values.  The  resulting  filter  is  Dgf.  Then, 
we  apply  a  filter  which  leaves  the  even  points  fixed  and  interpolates  to  get  the  odd 
points.  This  condition,  that  f  be  the  identity  on  even  points,  is  sufficiently  important 
to  warrant  a  separate  definition: 

Definition  2. 1 

The  lowpass  filter  f  is  said  to  be  an  'a  trous  filter  if  it  satisfies 

f2k  =  «(k)/V2  .  (2.9) 


The  V2  is  simply  a  convenient  means  of  including  the  normalization  factor  of  (1.6)  in 
the  filter.  The  result  of  the  entire  interpolation  operation,  as  pictured  in  figure  2.1,  is 
thus 


[  f*  .  (Dg')  ]„  =  (  F'Dg'  ]„ 

-SfJ_2kg(k)  (2.10) 

k 


k  k  ___  o  rj  #  t 

Noting  that  g(  — — n)  —  g( - )  and  inserting  the  approximation  (2.10)  into 


(1.6),  we  obtain 
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(2.11) 


w(2»  n)  *  E  fk— 2n-2m  gm  Sk 

k,ra 

=  £  ,  gn-m'  ^2m'  — k  sk 
k,  m 

=  [g  *  (A(f  *  s))]n 


which  is  simply  w'  of  equation  (1.7)  with  i  =  1.  Continuing  inductively  by  replacing  s 
in  (2.11)  with  s1_1,  we  find  that  w(2‘,2‘n)  ®  w„  for  all  i,  where  w„  is  given  by  (1.7), 
which  can  be  rewritten 


si+1  =  A  (f  *  s{)  (2.12a) 

w1  =  g  *  s1  .  (2.12b) 


Except  for  decimation  of  the  output  (the  undecimated  version  will  be  derived  in  sec¬ 
tion  3),  this  is  the  ' a  trous  algorithm  described  in  [4].  Thus,  we  see  that  the  'a  trous 
algorithm  is  simply  a  DWT  for  which  the  filter  f  (an  interpolator)  satisfies  condition 
(2.9)  and  the  filter  g  is  obtained  by  sampling  an  a  priori  wavelet  function  g(t). 

Remark  1 

The  definition  (1.6)  is  not  so  transparent  as  it  might  seem.  It  is,  of  course, 
intended  to  reflect  an  approximation  to  (1.5).  From  this  viewpoint  one  might  well 
consider  a  change  of  variables  t  —*  t/21  before  discretizing  (1.5).  Such  a  procedure 
certainly  alleviates  the  computational  problem  since  it  dilates  s(t)  (that  is,  samples  s  at 
21,  values  which  are  known)  rather  than  contracting  g(t).  However,  unless  the  original 
function  s(t)  was  highly  oversampled  (which  begs  the  computational  question),  the 
approximation  is  poor.  More  precisely,  to  accurately  approximate  s(t),  and  therefore 
also  (1.5),  we  must  sample  at  least  at  the  Nyquist  rate  rnyq  for  s.  Then  the  integral  for 
octave  i  requires  g(t)  to  be  sampled  at  a  rate  2lrnyq. 

Remark  2 

The  derivation  above,  as  well  as  that  in  [4],  of  the  'a  trous  algorithm  make  no 
statements  regarding  the  accuracy  of  the  approximation  (2.11)  or  even  of  the  discreti¬ 
zation  from  (1.5)  to  (1.6).  The  former  is  iterated  over  i  and,  hence,  to  succeed  must 
be  numerically  stable  in  some  sense.  This,  in  turn  depends  principally  on  choice  of 
the  filter  f.  A  major  step  towards  treating  this  question  lies  in  the  results  of  section  4, 
as  was  outlined  at  the  end  of  the  introduction.  Since  the  algorithm  is  exact  for  some 
g(t),  the  question  reduces  to  (a)  the  quality  of  the  approximation  g  ~  g  and  (b)  the 
effect  of  this  approximation  on  the  wavelet  integral  (1.1).  Inasmuch  as  g(n)  =  g(n) 
for  the  finite  set  of  integers  n  used  to  obtain  g  from  g(t),  it  is  sufficient  that  these 
functions  be  smooth  enough  and  decrease  fast  enough  at  infinity.  Conditions  on  f 
which  achieve  regularity  of  the  constructed  wavelet  g  are  found  in  [8].  Quantifying 
these  statements  remains  a  subject  for  future  study. 
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MULTIRESOLUTION  ALGORITHM 

Mallat’s  algorithm,  illustrated  in  figure  2.2,  has  essentially  the  the  same  tree 
structure  as  equations  (2.12)  (cf.,  [6]  —  [9]),  namely, 

si+1  =  A  (h  *  s’)  (2.13a) 

di+1  =  A  (g  *  s')  (2.13b) 

In  keeping  with  the  literature,  we  have  replaced  the  filter  f  with  the  filter  h,  which 
also  serves  to  indicate  that  this  class  of  filters  is  constrained,  as  detailed  below.  We 


Figure  2.2  The  Mallat  multiresolution  algorithm.  The  down- 
arrow  indicates  decimation. 


remark  that  none  of  these  filters  are  'a  trous  filters.  The  constraints  on  h  and  g  which 
insure  an  orthonormal  multiresolution  analysis  (cf.,  [6]  —  [9])  are 


?^2j— n  ^*2j - m  +  g2j  — n  §2j-m  ^nm  (2.14a) 

Sh2n_jfem_j  =  0  (2.14b) 

Egn  =  0  (2.14c) 

n 

£hn  =  V2  .  (2.14d) 

n 


Recalling  that  4  hj_j  and  that  Af  =  D,  we  may  rewrite  (2.14a)  and  (2.14b) 
as 

(HfD)  (AH)  +  (G+D)(AG)  =  I  (2.15) 

(AH)  (GfD)  =  0  .  (2.16) 
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Furthermore,  (2.15)  and  (2.16)  imply  (e.g.,  multiply  (2.15)  on  the  left  by  AH*) 


(AH)  (H*D)  =  I 

(AG)  (G*D)  =  I  (2.17) 


Thus,  H*D  and  G*D  are  injections  and  (2.13)  is  an  orthogonal  decomposition  of  the 
discrete  signal  s'.  That  is,  s1_1  =  HtDs'+GtDw1  with  the  scalar  product 
(HtDs‘)  •  (G*Dw‘)  =  0.  In  fact,  equation  (2.15)  is  a  paradigm  for  inverting  the 
transform.  These  concepts  are  illustrated  in  figure  2.3. 


..i-i 


AH 


s1 


AG 


d* 


* 


G*D 

d’ 


Figure  2.3  Illustration  of  a  single  stage  and  its  inverse  in  the 
multiresolution  algorithm  for  orthonormal  wavelets. 


Furthermore,  from  (2.14)  it  follows  that  (2.13)  represents  a  wavelet  decomposi¬ 
tion  (see  [6],  [8],  [9])  as  follows:  There  exists  a  scaling  function  <£( t)  whose  Fourier 
transform  is  given  by 


<£(w) 


1 

V2 


hz(4) 

zv  2, 


(2-18) 


Expression  (2.18)  forms  the  basis  of  a  recursion, 


(2.19) 


which  in  the  time  domain  takes  the  form 


0(t)  =  E  h_k  V2  0(2t  -  k)  .  (2.20) 

k 

With  some  additional  computation  (cf.,  section  4),  one  may  demonstrate  that  the 
translates  and  dilates  of  <f>. 
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(2.21) 


*«  -  w  “ n) 

have  the  property 

4+I(0  =  S[AH]nk^(t)  .  (2.22) 

Note  that  the  above  definitions  differ  in  the  sign  of  i  from  those  of  [8]. 

Finally,  define 

m  &  =  s  V2  g_k  —  k)  .  (2.23) 

Then,  using  (2.14)  and  the  above  properties  of  (f>,  one  can  show  that  the  family  of 
wavelets/ 

«(•)  4  *<7  -  ■>)  (2-24) 

are  orthonormal  (  /  t/>„(t)  4(0  dt  =  *j4.k  ),  and  that  the  d1  are  the  coefficients  of  the 
expansion  of  s(t)  in  terms  of  the  ip'n. 

More  precisely,  the  integer  translates  <j>(i  —  k)  of  the  scaling  function  form  a 
basis  for  L2(R),  and 

4  =  s(t)  Mjr  -n)  =  f  s(t)  4(0  (2.25) 

provided 

Sk  =  /  s(t)  ^(t-k)  .  (2.26) 

Then,  if  s(t)  is  in  L2(R),  the  expansion  of  s(t)  is 

s(t)  =  ?  d*  4(0  •  (2-27) 

i  n 

This  follows  from  (2.25)  and  orthonormality  since  completeness  of  the  <t>( t  —  k) 
implies  completeness  of  the  0„(t)  . 

1  In  contrast,  the  gj,(t)  -  ■  g(-^--n)  of  (1.5)  are  not  generally  orthogonal.  It  is  the  filter 

constraints  (2.14a)  and  (2.14b)  that  insure  orthonormality.  Dropping  these  two  constraints  in 
section  4,  we  develop  a  structure  identical  to  (2.21)  —  (2.24);  however,  the  constructed  gj,(t) 
need  not  be  orthogonal. 
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Until  recently,  the  only  known  wavelets  with  compact  support  (i.e.,  zero  outside 
a  finite  interval)  were  the  Haar  functions,  generated  by  —  1  for  0  t  <  1/2;  -1 
for  —1/2  <  t  <  0;  and  0  elsewhere.  Ingrid  Daubechies  (cf.,  [8],  [19])  has  uncovered 
an  entire  family  of  finite  length  filters  satisfying  (2.14),  demonstrating  that  the 
corresponding  wavelets  defined  by  (2.18),  (2.23)  and  (2.24)  are  orthonormal  as  a 
consequence  of  (2.15)  -  (2.16)  and  have  compact  support  since  they  were  generated  by 
finite  length  filters  (cf.,  section  4,  equation  (4.10))  The  first  two  of  these  filters  are 

h+  =  ^2  (1’  V  ’  (2.28a) 

and 

hf  =  (1+V3,  3+V3,  3-V3,  1-V3)  (2.28b) 

where  the  first  component  of  h  is  on  the  left.  The  wavelets  corresponding  to  (2.28a) 
are  exactly  the  Harr  function  mentioned  above. 

Some  additional  remarks  relating  the  two  algorithms  are  in  order.  The  condi¬ 
tions  (2.14c)  and  (2.14d)  effectively  make  g  a  bandpass  filter  and  h  a  lowpass  filter 
(e.g.,  an  interpolation  filter)  with  the  sum  on  gk  analogous  to  the  condition 
/g(t)dt  =0.  Also,  d1+1  and  w1  correspond;  the  additional  decimation  appearing  in 
(2.13b)  would  not  appear  in  a  translation  invariant  version  of  Mallat’s  algorithm  (cf., 
section  3).  On  the  other  hand,  although  the  discrete  filters  g  play  algorithmically 
identical  roles,  the  g„(t)  of  the  'a  trous  algorithm  are  not  the  wavelet  vectors  of  a  func¬ 
tional  expansion.  Rather  the  g^(t)  and  0n(O  are  the  duals  of  a  set  of  vectors  for 
which  coefficients  of  the  signal  expansion  are  W2n  and  d„  respectively  (cf.,  [3],  [15]). 
That  is,  they  are  the  coefficients  of  an  expansion  of  the  form  s(t)  =  E 
<s,  $„(t)  where  <  >  indicates  the  L2  inner  product,  and  is  the  dual  basis  oi,n 

frame  (see  [3])  of  ip'n.  In  Mallat’s  algorithm,  since  the  Vi(t)  are  orthonormal,  the 
basis  and  its  dual  coincide.  Thus,  in  many  senses,  the  discrete  filters  g  are  more  fun¬ 
damental  than  the  wavelets  themselves.  Even  in  the  orthonormal  case,  it  is  usually 
the  coefficients  which  are  of  major  interest;  the  actual  wavelets  0„(t)  or  the  duals  of 
the  g„(t)  are  rarely  computed. 

Finally,  in  anticipation  of  section  5,  let  us  form  the  squares  of  the  two  previous 
examples 


ht  *  h  =  -yjd,  1)  *  1)  =  (1/2,  1,  1/2)  (2.29a) 

and 

h’*h  (229b) 

where  *  indicates  convolution,  and  h  given  by  (2.28a)  and  (2.28b)  respectively.  Note 
that  the  filters  f  =  h+  *  h  /V2  are  a  trous  filters  and  that  the  interpolation  equation 
(2.10)  holds  exactly  if  g(t)  is  a  polynomial  in  t  of  degree  one  or  three  respectively. 
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3.  WAVELET  TRANSFORMS  WITHOUT  DECIMATION 


As  has  often  been  pointed  out  [7],  the  recursions  (2.12)  and  (2.13)  are  not,  in 
general,  translation  invariant.  In  contrast,  the  original  undecimated  'a  trous  algorithm 
(cf.,  [4],  [5]),  which  is  pictured  in  figure  3.1  and  consists  entirely  of  convolutions,  is 
patently  translation  invariant.  In  this  subsection  we  shall  use  that  property  to  provide 
a  formal  definition  for  an  undecimated  discrete  wavelet  transform  w,  then  demon¬ 
strate  that  w^,  =  w}^,  and  also  show  that  w  is  computed  by  the  algorithm  of 
figure  3.1. 


D'f 


pi-f  1 


D‘+1f 


,i+2 


D'g 


D;+1g 


w 


w 


i+1 


Figure  3.1  The  (undecimated)  discrete  wavelet  transform. 
The  filters  D’f  are  obtained  from  f  by  inserting  2*— 1  zeros 
between  each  pair  of  filter  coefficients.  The  operation  of 
filtering  is  understood  to  mean  convolution. 


Let  Tra  be  the  operation  of  translation  by  m;  i.e., 

(Tms)n  A  Sn_m  (3.1) 

In  order  to  include  the  dependency  of  w1  on  s°,  we  add  it  as  an  argument,  writing 
w'(s°).  Equations  (2.12a)  and  (2.12b)  become 

w'(s°)  =  G(AF)‘  s°  .  (3.2) 

(Recall  that  Gjj  =  gi-j-)  Finally,  we  shall  have  need  of  the  following  important  identi¬ 
ties,  which  are  proved  in  appendix  A:  For  any  F  of  the  form  Fjj  =  fj_j, 

Lemma  3.1 


and 


(AF)^  =  (AF)'0  k_2i, 


(3.3) 
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Lemma  3.2 


i—i 

S  (AF)‘nk  eikw  =  ei2inw  11  tz(2'u)  .  (3.4) 

k  j— 0 

As  expected,  w1  is  not  translation  invariant, 

[»‘(T»s0))„  =  E  [G(AF)‘|„k  4_m 
=  E  [G(  AF)']nk+m  sj 

#  E  (G(AFy]0_„,k  s?  .  (3.5) 

For  example,  [AF];j  =  f2i-j  /  [AF]0  ;_j.  However,  if  we  replace  m  by  2‘m  in  (3.5) 
and  use  Lemma  3.1,  the  last  step  becomes  an  equality  so  that 

[»ia'*ms0)]„  =  [*i(s°)l„-„  .  (3.6) 


Thus,  translating  s°  by  2‘m  translates  octave  i  by  m. 

Note  that  the  zeroth  element  of  a  series  is  invariant  under  decimation  so  that  w„ 
and  wj,  should  coincide  at  n  =  0.  Using  this  fact,  we  obtain  the  nth  output  of  the 
undecimated  discrete  wavelet  transform  by  translating  the  signal  back  by  n  samples 
and  taking  the  decimated  transform  at  time  zero.  More  precisely, 

Definition  3.1 

Define  the  undecimated  discrete  wavelet  transform  w  in  terms  of  the  decimated 
transform  w  by 


wi  4  [^(s0)],,  A  [wi(T_ns°)]0  ,  (3-7) 

We  see  that  the  desired  invariance  is  achieved, 

[wi(Tms°)]n  =  [w'(T_nTms°)]0 

=  [wi(Tni_s0)]o 

=  [«*(•*»)]„  -  (3.8) 

It  is  also  clear  from  (3.6)  and  (3.7)  that  sampling  wj  every  2'  points  produces  exactly 
w„,  that  is, 


w; 


w 


2"n 


(3.9) 
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Next,  we  show  that  w  may  be  computed  by  the  filter  sequence  pictured  in 
figure  3.1.  The  proof  is  obtained  by  taking  z-transforms.  From  (3.2)  and  (3.7) 


Ewje-i""  =  E  [G(AF)i]0„  s„„  e-‘ 

n  mn 


=  £  [G(AF)']0m  eim"  sz(cj)  . 

m 


(3.10) 


Applying  Lemma  3.2  to  (3.10),  we  have 


wz(w)  4  £  w^,  e  ir 


i— 1 


=  gz(2^)  I!  fz(2iw)  sz(u) 
j-o 


(3.11) 


where  i  =  0  is  understood  to  mean  there  are  no  factors  of  fz.  As  described  in  the 
next  paragraph,  this  is  exactly  the  z-transform  of  the  algorithm  pictured  in  figure  3.1. 

It  is  easy  to  see  that  D‘f  is  f  with  21— 1  zeros  inserted  between  every  pair  of  filter 
coefficients  and  that  its  z-transform  is  fz(2‘o;).  That  is, 

n  =  2'm 

0A  other  <312> 

and 

(D‘f)z  =  (z(Tu)  .  (3.13) 

Equation  (3.11)  is  then  equivalent  to 

si+1  =  (D’f)  *  s*  (3.14a) 

w1  =  (D‘g)  *  s'  (3.14b) 


where  s°  4  s.  This  is  essentially  the  original  (undecimated)  'a  trous  algorithm  found 
in  [4]  and  [5].  However,  we  would  like  to  emphasize  that,  since  the  development  in 
this  subsection  has  not  made  use  of  any  filter  constraints,  the  general  equivalence  of 
the  decimated  output  of  the  algorithm  pictured  in  figure  3.1  (equations  (3.14))  and 
that  of  figure  1.1  (equations  (2.12))  follows  for  arbitrary  filters. 

An  alternative  implementation 

A  second  possibility  for  the  implementation  of  w  is  to  use  the  algorithm  in  figure 
1.1  and  proceed  directly  from  equation  (3.7).  That  is,  to  output  w'n,  one  simply 
translates  the  signal  by  n  and  then  computes  v/'0.  As  an  algorithm,  this  takes  the  form 
(a)  compute  w),;  (b)  translate  the  signal  by  one  sample;  (c)  go  to  step  (a).  Moreover, 
in  view  of  (3.6),  one  need  not  reperform  the  entire  recursion  (i.e.,  equations  (2.12)) 
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for  every  time  point  n  in  order  to  obtain  w'.  Rather,  at  each  octave,  the  decimation  is 
replaced  by  a  split  into  even  and  odd  sequences,  each  of  which  is  a  starting  point  for 
the  next  octave  (cf.,  figure  3.2).  A  few  examples  suffice  to  convince  one  that  if,  at 
octave  i,  n  mod  2‘  =0,  then  one  takes  the  upper  branch;  if  n  mod  2‘  =  1,  then  one 
takes  the  lower  branch.  A  rigorous  derivation  follows  from  the  formula  (cf.,  [20]) . 


ATms 


m 

T  2  seven  m  even 

m— 1 

T  2  sodd  m  odd 


(3.15' 


(w1),,, 


Figure  3.2  Diagram  of  an  implementation  of  the  undecimated  DWT. 


We  remark  that  figures  5  and  6  are  computationally  equivalent  provided  that  the 
algorithm  in  figure  3.1  is  implemented  efficiently.  The  code  must  be  written  so  as  to 
omit  multiplication  by  the  zero  elements  of  filters  D‘f.  (They  are  mostly  zeros  for  i  > 
2.)  Depending  on  the  number  of  octaves,  the  computational  burden  still  remains 
much  greater  than  that  of  the  decimated  algorithm  (i.e.,  figure  1.1);  however,  there  is 
considerable  parallelism  which  may  be  sufficiently  exploited  on  suitable  hardware  to 
produce  a  realtime  implementation.  [5] 
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4.  THE  DWT  AS  AN  EXACT  WAVELET  TRANSFORM 


Regardless  of  the  filters  used,  one  can,  of  course,  perform  the  recursions  (2.12) 
or  (2.13)  on  the  sampled  signal  s.  Moreover,  provided  that  f  (respectively,  h)  is 
lowpass  and  g  bandpass,  the  procedure  may  be  interpreted  physically  as  a  bank  of 
proportional  bandwidth  filters  (cf.,  [21]  —  [24],  also  section  6).  In  the  present  section, 
we  examine  the  mathematical  significance  of  relaxing  the  filter  constraints  (2.9)  and 
(2.14).  Our  goal  will  be  to  relate  the  more  general  filter  bank  to  the  continuous 
wavelet  transform,  thus,  in  a  sense,  justifying  the  term  DWT.  In  this  endeavor,  the 
major  questions  we  shall  address  are:  for  what  functions  g(t)  is  the  recursion  (2.12)  an 
exact  implementation  of  (1.6)  and  for  which  g(t)  and  s(t)  do  (1.5)  and  (1.6)  coincide? 
The  general  answer  is  that  we  are  able  to  construct  such  a  g  provided  the  discretized 
signal  lies  in  the  appropriate  subspace  (cf.,  equation  (2.26)).  A  somewhat  surprising 
result  is  that  it  is  necessary  and  sufficient  for  f  to  be  'a  trous  for  condition  (2.26)  to  be 
dropped.  Our  approach  shall  be  to  mimic  the  construction  of  orthonormal  wavelets 
outlined  in  section  2. 

CONSTRUCTION  OF  THE  SCALING  FUNCTION  <t> 

We  begin  with  the  stipulation  of  the  existence  of  a  scaling  function  <f>(t )  with 
Fourier  transform 


0(w)  4  fl 

j-i 


1  fz(^) 

zv  2,  > 


V2 


(4.1) 


where  f z(w)  =  (ff)z(w)  is  the  z-transform  of  ff.  To  emphasize  the  nonorthogonality 
of  the  corresponding  wavelets,  we  retain  the  symbol  f  rather  than  h.  Note  that  this 
function  6  need  not  have  (and  in  general  does  not  have)  all  of  the  properties  of  the 
orthonormal  6  outlined  at  the  end  of  section  2. 

For  (4.1)  to  converge  to  a  nonzero  function,  the  factors  must  approach  one. 
Thus,  fz(0)  =  1,  which  implies 


E  fk  =  V2  .  (4.2) 

Even  though  6  could  be  normalized  differently  by  the  inclusion  of  a  factor  in  (4.1), 
the  filter  f  must  necessarily  obey  the  lowpass  condition  (4.2);  i.e.,  (2.14d).  Note  also 
that,  under  the  chosen  normalization,  /  <^(t)dt  =  6(0)  =  1.  However,  without  some 
additional  conditions,  the  relationship  of  <£(t)  to  6(^)  remains  somewhat  tenuous. 
Even  under  pointwise  convergence,  the  limit  may  be  a,  highly  discontinuous,  fractal 
function.  [9]  Suitable  regularity  conditions  for  the  inverse  Fourier  transform  of  a  pro¬ 
duct  of  the  form  (4.1)  to  converge  to  a  reasonably  behaved  (e.g,  L^R),  L2(R),  and/or 
continuous)  function  may  be  found  in  [8]  and  [25].  The  results  are  summarized  in 
appendix  B. 

A  time  domain  representation  of  6(t)  in  terms  of  f  can  be  derived  from  (4.1). 
Let  \  be  the  indicator  function  of  [—1/2, 1/2) 
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fl  for  t  6  [-1/2,  1/2) 
X(0  -  jo  other 


(4.3) 


From  (4.1)  and  Lemma  3.2,  it  is  not  difficult  to  see  that 

=  lim  £  ((AF)i)0k  V2iX(2h-k)  •  (4.4) 

j— *oo  k 

g 

In  fact,  the  Fourier  transform  of  (4.4)  is 


lim  £  (AF)ok  e 

i— >oc  k 


1  U  , 

= 


lim  FI 

I— »O0  j«T 


which  is  just  (4.1).  Note  that  the  existence  of  the  function  (4.4)  could  have  been 
taken  as  the  starting  point  for  our  analysis,  since  our  proofs  will  not  make  use  of 
(4.1). 

We  continue  our  parallel  with  the  development  of  Mallat’s  algorithm  in  section  2. 
Define  <f>'n(t)  by 

Definition  4.1 

K(0  -  -ft  -  ")  •  («) 

Then,  substitution  of  2't  —  n  for  t  in  (4.4)  and  the  use  of  Lemma  3.1  yield 

V2j  4>{2'l  -n)  =  lim  £  ((AF)Ook  V2'+i<£(2i+it  -  2'n  -  k) 

j— «x>  k 

=  Hm  £  ((AF)H)  V2^(2H-k) 

j— *oo  k  ’ 

=  lim  £  ((AF)i-i )nk  V2iX(2H-k)  .  (4.7) 

j— *oo  k 

On  replacing  i  by  — i,  this  becomes 

0*(t)  =  lim  E  ((AF)i+i)nk  V^t-k)  .  (4.8) 


8 


This  equality  is  immediate  from  Lemma  3.2  and  y(w) 
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An  immediate  consequence  of  (4.8)  is 

C1  -  ?  IAF]nk  «  (4.9) 

k 

paralleling  (2.22).  Thus,  we  see  that,  despite  their  lack  of  orthogonality,  the  $,(t) 
have  retained  most  of  their  structure.  Furthermore,  if  f  is  a  finite  filter,  then  <?i(t)  has 
finite  support.  [8]  More,  precisely,  suppose  the  coefficients  of  f  are  zero  outside 
[— N_,  N+].  Let  7?J( t)  4  ((AF)j)0k  V2j  x(2Jt  —  k).  Then,  t)  converges  to  <£(t),  and, 
as  in  (4.7),  we  have 


?/j+1(t)  =  S  fkf  ??J(2t  -  k)  .  (4.10) 

k 

We  find,  for  example,  that  rj1  is  zero  for  t  <  tj  where  tj  =  (tj—1  —  N_)/2.  With 

t0  =  — 1/2,  and  with  a  similar  calculation  for  the  right  half  interval,  it  follows  that 

<£(t)  =  lim  ^(t)  is  zero  outside  [— N_,  N+). 
j— 00 

EXACTNESS 

To  avoid  confusion  and  stress  their  differences,  let  us  first  recapitulate  some 
definitions.  Four  different  transforms  W(a,  b),  W(2‘,  2'n),  w(2',  2'n),  and  w^  have 
been  mentioned  ((1.1),  (1.5),  (1.6)  and  (2.12)).  We  retain  a  terminology  parallel  to 
Fourier  transforms,  namely,  wavelet  transform  (WT),  wavelet  series  ,  discretized 
wavelet  series,  and  discrete  wavelet  transform  (DWT)  respectively.  The  first  two 
transforms  involve  integrals  of  a  continuous  signal;  the  latter  two  contain  sums  of 
sampled  signals.  The  first  three  utilize  a  continuous  wavelet  function  g(t),  the  last  one 
employs  the  discrete  filters  g  and  f.  For  consistency,  we  shall  continue  our  develop¬ 
ment  using  decimated  transforms;  however,  the  results  hold  without  change  for  unde¬ 
cimated  transforms.  This  follows  immediately,  since  they  coincide  for  n  =  0,  and  the 
undecimated  transforms  may  be  obtained  at  time  n  =  n0,  by  translating  the  signal  by 
n0  samples  and  taking  the  transform  at  n  =  0.  (See  section  3  where  definition  3.1 
remains  valid  for  W(2‘,  n)  and  w(2‘,  n)). 

Our  starting  point  shall  be  a  signal  s(t)  and  discrete  filters  f  and  g  with  w1  defined 
by  (2.12)  or  d'  by  (2.13),  i.e., 


si+1  =  AF  s! 

w!  =  G  s'  (4.11) 

di+1  4  Aw'  . 

Recall  that  the  matrices  Fjj  and  Gjj  are  given  by  f;_j  and  g;_j  respectively.  Of  course, 
we  must  also  specify  an  initialization  of  the  recursion  (4.11)  for  some  i;  for  example, 

^  At  times,  we  shall  prefer  the  term  sampled  WT  rather  than  wavelet  series  in  order  to  em¬ 
phasize  its  role  as  a  restriction  of  the  continuous  transform. 
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for  the  zeroth  octave  s°.  The  obvious  choice  is 


s°  £  s(n)  ;  (4.12a) 

however,  we  shall  also  consider 

s°  4  E  d(k  —  n)  s(k)  (4.12b) 

k 

which  relates  to  the  discretized  wavelet  series  w(2‘,  2’n),  and 

s „  =  f  4*  —  n)  s(t)  dt  (4.12c) 

which  corresponds  to  the  sampled  WT  (wavelet  series).  For  a  given  g,  we  shall  con¬ 
struct  a  continuous  function  g(t)  such  that  the  DWT  of  equation  (4.11)  is  an  exact 
implementation  of  the  discretized  wavelet  series  under  (4.12b)  and  of  the  wavelet 
transform  under  (4.12c). 

Define  g(t)  by 

g(t)  ^  E  <t>(t  +  k)  gk  =  E  d>(l  -  k)  gl 

and  gj,(t)  by 

?*(•;  =  ^rg(ir  “n)  • 

Then,  using  (4.6),  (4.9),  (4.13),  and  (4.14),  we  have 

gn(0  =  ^7-  g(^T  ~n) 

=  E  g l  —y-r  <t>{—  —  n  —  k) 
k  V21  2‘ 

=  ee  gn_k  m 

k 

=  E  [G(AF)‘]nk  4°(t)  .  (4.15) 

k 

If  we  take  (4.12b)  as  the  definition  of  s°,  the  discretized  wavelet  series  (equation 
(1.6))  takes  the  form 

w(2\  2'n)  4  E  g„(m)  s(m) 

mk 

=  [EG  (AF)']nk  4°(m)  s(m) 

m 

00 


(4.13) 

(4.14) 


-  [G  (AF)1  s°]„ 


(-1.16) 


which,  by  (4.11)  is  exactly  w„  .  Furthermore,  under  (4.12c)  we  have  (again  using 
(4.15)) 

W(2‘,  2'n)  &  {  g‘(t)  s(t)  dt 

=  /  E[G(AFy]nk40(t)s(t)dt 

k 

=  E[G  (AF)’)nk  /  <^k(t)  s(t)  dt 

k 

=  [G  (AF)'  s°]n  ,  (4.17) 


again  w'n  . 

Finally,  let  us  investigate  the  significance  of  the  'a  trous  condition;  i.e.,  of  the 
constraint  (2.9),  ^2k  =  (5k0/V2.  We  prove  the  following  theorem: 

Theorem  4.1 

f  is  an  'a  trous  filter  <*>  4>(n )  =  <5n0. 

Proof 

Letting  i  =  —1,  n  =  0,  and  t  =  n  in  equation  (4.9)  gives 

H n)  =  E  (AF)0k  V2  <j>(2n  -  k)  .  (4.18) 

k 

Then, 

<Kn)  =  <5n0  =>  V2  f2n  =  6n0  .  (4.19) 

Conversely,  suppose  that  f2n  =  <5n0/V2.  Then,  since  x(~k)  —  ^ko*  equation  (4.8) 
with  i  =  0  and  t  =  0  implies 


<£(— n)  =  lim  (AF)in0  Vl> 

J-.0O 

=  lim  E(AF)i-1nk(AF)k0V2V2i-1 

j— *oo  k 

=  lim  (AFy_1n0  V2i-‘ 

J—OO 

.  .  .  =  lim  (AF)n0V2  =  6n0  .  (4.20) 

J-*00 

The  import  of  theorem  4.1  quickly  follows.  Equation  (4.12b)  implies  (4.12a)  for 
arbitrary  signals  if  and  only  if  <£(n)  =  <5n0.  Thus, 
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Corollary  4.1 

The  algorithm  (4.11)  is  an  exact  implementation  of  the  DWT  with  s°  set  equal  to 
the  sampled  signal  if  and  only  if  f  is  an  'a  trous  filter. 

Furthermore,  g  is  also  related  to  g(t)  by  sampling: 

Corollary  4.2 

If  f  is  an  'a  trous  filter,  then  g(t)  defined  by  (4.13)  satisfies  g(n)  =  g„. 

It  remains  to  investigate  the  role  of  dl.  Define 

V*(t)  4  V2  g(2t)  =  V2Sgk  <K2t  +  k)  (4.21) 

k 

which  coincides  with  (2.23).  Then,  provided  s°  satisfies  (4.12c),  (4.11),  and  (4.17) 
imply 

d*  =  (AsS-% 

=  win1 

=  /  SnHO  S(t)  dt 

-  /  ^  2  }  S(t)  dt 

=  /  mi,  (t)  s(t)  dt  .  (4.22) 

Hence,  although  the  (i/^)h(t)  are  not  orthogonal,  the  Mallat  algorithm  still  computes 
the  wavelet  transform.  Of  course,  under  (4.12a)  and  'a  trous,  or  for  (4.12b),  we  have 
the  counterpart  of  (4.16) 

d'  =  E  (#)j,(k)  s(k)  .  (4.23) 

It  is  interesting  that,  in  a  sense,  the  decimated  wavelet  transforms  (1.5b)  and  (1.6) 
contain  superfluous  information.  That  is,  they  are  underdecimated  by  a  factor  of 
two,  and,  thus,  w1  properly  belongs  to  octave  i  +  1. 

SUMMARY 

Let  us  summarize  the  results  of  this  section.  We  are  given  discrete  filters  f  and  g 
such  that  (4.4)  is  well  defined.  Define  g(t)  by 

g(t)  =  £  <f>{\  -  k)  (4.24) 

k 
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with  corresponding  transforms  sampled  WT  (wavelet  series) 

W(2\  2‘n)  =  /  g‘(t)  s(t)  dt 


(4.25) 


and  discretized  wavelet  series 

w(2‘,  2*n)  4  E  g‘(k)  s(k)  .  (4.26) 

Let  <f>n  •  s  stand  for  the  scalar  product  E  <^n(k)sk,  and  ^n(t)  •  s(t)  for  the  L2  scalar  pro¬ 
duct  /  <An(t)  s(t).  Then: 


f  is  'atrous  =>  g(n)  =  gn 

(4.27) 

For  s  discrete: 

Sn  =  4>n  "  S  => 

w(2‘,  2‘n)  =  w‘ 

(4.28a) 

s°  =  s  and  f  is  ' a  trous  => 

w(2’s  2’n)  =  wj,  . 

(4.28b) 

For  s(t)  continuous: 

Sn  =  </>n(0  •  S(t)  =£> 

W(2\  2‘n)  =  wj,  . 

(4.29) 

Define  tf*(t)  4  V2  g(2t).  Then  the 
replaced  by  #(t);  i.e., 

corresponding  properties  hold  for  d’ 

with  g(t) 

Sn  =  ‘  s(0  => 

<  =  f  m  (t)  s(t)  . 

(4.30) 

S„  =  (j)n  •  s  or 

s°  =  s  and  f  is  'a  trous  => 

di  =  E  (^)i(k)  s(k)  . 

k 

(4.31) 

The  above  results  extend  immediately  to  the  undecimated  transforms, 
w(2‘,  n),  and  wj,  by  translation  invariance  and  Definition  3.1. 

W(2‘,  n), 
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5.  LAGRANGE  INTERPOLATION  FILTERS 


A  reasonable  class  of  'a  trous  interpolators  to  consider  for  equation  (2.10)  are 
those  that  are  exact  for  polynomials  P(t)  of  degree  <  M  for  some  M,  i.e.,  for  which 

-  S  fJ_,„P(k)  .  (5.1) 

For  reasons  that  will  become  clear  very  shortly,  we  shall  call  these  filters  Lagrange  'a 
trous  filters.  Since  the  'a  trous  filter  f  satisfies  f2ic  =  <50k/V2,  equation  (5.1)  is  an 
identity  for  n  even.  Let  a  contain  the  odd  components  of  ff, 


3k  = 


V2  f^k-i  for  k  >  0 
V2  f}k+1  for  k  <  0  . 

1  for  k  =  0 


Then  (5.1)  is  equivalent  to 


P(f)  =  fn  ak  P(~—k)  +  E  akP(ili-k)  nodd 

L  k>0  l  k<0  l 


(5.2) 


(5.3) 


for  all  polynomials  P  of  degree  <  M.  We  proceed  to  express  the  ak  in  terms  of 
Lagrange  polynomials,  and  to  show  that  the  above  conditions  arc  essentially 
equivalent  to  f  =  h  *  hf  where  h  is  an  appropriate  Daubechies  filter. 

CONSTRUCTION  OF  THE  FILTER  a 

First,  we  parameterize  the  family  of  filters  a  satisfying  (5.3)  by  M  +  l,  the  dimen¬ 
sion  of  the  space  of  polynomials  for  which  it  must  hold.  For  such  a  relationship  to 
exist,  one  must  relate  the  length  of  the  filter  (the  number  of  unknowns)  to  M.  To 
accomplish  this,  we  shall  assume  that  a  has  exactly  the  minimum  number  of  coeffi¬ 
cients  needed  to  satisfy  equation  (5.3)  .  We  further  assume  that  a  has  symmetric  sup¬ 
port;  i.e.,  there  is  an  N  such  that  ak  =  0  for  |k|  >  N  and  ak  ^  0  for  |k|  =N.  This 
assumption  is  not  unreasonable,  at  least  for  symmetric  wavelets  g(t),  since  there  is  no 
a  priori  reason  to  distinguish  between  t  and  -t,  and  one  would  even  expect  a  to  be 
symmetric.  We  shall  see,  in  fact,  that  the  weaker  condition  of  symmetric  support 
joined  with  the  previous  constraints  implies  that  a  actually  is  symmetric. 

The  minimum  number  of  coefficients  (unknowns)  equals  the  dimension  of  the 
space,  i.e.,  satisfies  2N  =  M  +  l.  Thus,  the  sums  in  (5.3)  go  from  1  to  N.  More¬ 
over,  since  for  any  n,  the  polynomials  Q(x)  4  P(x  —  (n— 1)/2)  also  form  a  basis, 
equation  (5.3)  is  equivalent  to 

P(^-)  =  £  akP(l— k)  +  E  akP(— k)  nodd  ,  (5.4) 

2  k  =  l  k=«-l 
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which  must  hold  for  all  P(x)  of  degree  <  2N— 1. 

We  pick  out  the  k*h  coefficient  by  letting  P  be  the  Lagrange  polynomial 


n  (x  -  i) 

Lj2N_1(x)  4  t y  >,  j  in  [-N+1,  N]  .  (5.5) 

Mi 

Then,  L2N-1(k)  =  <5jk,  so  that  replacing  P  in  equation  (5.4)  with  L2N_1,  we  get 

aj  =  L^-^y)  for  j  =  1,  ...,  N  (5.6a) 

a_j  =  Lj2N_1(y)  for  j  =  1,  ...,  N  .  (5.6b) 

Inasmuch  as  the  Lagrange  polynomials  L2N_1  form  a  basis  for  polynomials  of  degree 
less  than  or  equal  to  2N  —  1,  these  ak  are  in  fact  the  unique  solution  to  (5.4). 

It  is  also  straightforward  to  see  from  (5.5)  and  (5.6)  that  a  is  symmetric,  i.e.,  for 

j>0 


n  (i/2 -i) 

=  _ 

n.(-j+i-i) 

n*i-j 


n  (-1/2  i) 
_ 


n  (-j  +  i) 


n  (i/2  -  i) 

m _ 

n  (i  -  i) 

>*j 


(5.7) 


In  summary,  we  have 


Theorem  5.1 


Let  f  be  an  a  trous  filter,  i.e., 


^2k  — 


V2 


'Ok 


(5.8) 


Assume,  furthermore,  that  f  is  real  with  symmetric  support  described  by 
k  e  (— 2N+1,  2N  — 1].  Then,  f  is  a  Lagrange  'a  trous  filter,  that  is,  (5.1)  holds  for  all 
polynomials  P  of  degree  <2N  — 1,  if  and  only  if  the  odd  components  of  f  are  deter¬ 
mined  by  equations  (5.6).  Furthermore  f  is  necessarily  symmetric. 


RELATIONSHIP  TO  DAUBECHIES  (QMF)  FILTERS 

In  [8],  Ingrid  Daubechies  constructs  essentially  the  entire  class  of  finite  length 
filters  h  which  satisfy  (2.14)  and  fulfill  suitable  regularity  conditions  on  (2.18).  Expli¬ 
citly,  they  take  the  form 

h(z)  =  (yd  +*))NQ(z)  (5-9) 


where  h(z)  4  Ehnz  n  is  the  z  transform  of  h,  and  Q  is  an  appropriately  constrained 


27 


polynomial.  (In  this  section  it  is  convenient  to  express  the  z  transform  as  a  polyno¬ 
mial  which  we  denote  h(z)  where  hz(w)  =  h(e,w).)  The  derivation  uses  a  specific  g, 
which  up  to  a  phase  factor  is  given  by 

gn  =  (-l)"h(l  -n)  ,  (5.10) 

and,  as  a  consequence,  (2.14a)  reduces  to 

Ah  hf  D  =  I  .  (5.11) 

In  other  words,  [hh^n  =  S0n,  which  is  the  'a  trous  condition  for  hl^/VE  Finally,  if 
Q  is  taken  to  be  of  minimal  degree  (which  turns  out  ro  be  N  —  1),  | O | 2  is  unique.  [8] 
In  other  words,  the  squares  of  these  filters  are  characterized  completely  by  satisfying 
(2.14d)  and  (5.11),  and  being  of  the  form  (5.9)  with  degree  2N  —  1.  We  proceed  to 
show  that  the  h  *  hT  equal  the  Lagrange  'a  trous  filters. 

Let  h  be  the  Daubechies  filter  of  order  2N.  Since  h  *  h*  is  symmetric,  the  above 
conditions  are  equivalent  to 


N  N 

B(z)  £  h  *  hf  =1  +  E  bk  z2k-x  +  Ebk  z-2k+1  (5.12a) 

and 

B(z)  =  y(y(l+z))N({(l+z-1))NQ(z)Q(z-1)  .  (5.12b) 

That  is,  V2  f  =  h  *  h+  if  and  only  if  f(z)  is  of  the  form  (5.12a)  -  (5.12b)  and  is  of 
degree  2N  —  1  in  z  (respectively,  in  z-1). 

Next,  we  show  that  the  bk  coincide  with  the  ak  of  equation  (5.3).  We  multiply 
(5.12a)  by  zn  for  an  arbitrary  integer  n, 


zn  B(z) 


+ 


N 

V 


bt  z 


n-l+2k 


+ 


N 

E 

i 


,n+l-2k 


The  2N  zeros  at  z  =  -1  in  (5.12b)  imply  that 


/nznB(z)) 


and  since  B(— 1) 


0,  we  also  have 


i  =  1,  ...,  2N  — 1 


(5.13) 


(5.14a) 
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1 


(5.14b) 


N 

2  E  bk  = 

l 


Next,  we  note  that 

=  m  (m  —  1)  •  •  •  (m  —  i  +  1)  (— l)m_1  for  i  >  0.  (5.15) 

Define  a  set  of  polynomials  P;  by 

P;(x)  4  2x  (2x  —  !)•••  (2x  —  i  + 1)  .  (5.16) 


ffl(zm) 

(dzy 


Finally,  setting  the  ith  derivative  of  the  right  hand  side  of  equation  (5.13)  at  z  =  —  1  to 
zero,  using  (5.15)  and  definition  (5.16),  we  obtain 

Pi(f)  =  E  bkPi(^-+k)  +  ?  bkPj(Jl±l-k)  (5.17) 

Z  k>0  Z  k>0  Z 

where  i  =  1,  ...,  2N— 1.  The  signs  work  out  since  n— l+2k— i  and  n+l-f2k  —  i  have 
the  same  parity  while  n— i  differs.  Define  Po(x)  =  1.  Then,  since  the  polynomials 
P0(x)  and  P;(x)  for  i  =  1,...,  2N— 1  form  a  basis  for  polynomials  of  degree  <  2N— 1, 
and  since  equation  (5.14b)  implies  (5.17)  for  i  =  0,  equation  (5.17)  must  hold  for  arbi¬ 
trary  polynomials  P(x)  of  degree  <  2N  —  1.  Replacing  P;  by  P  in  equation  (5.17)  and 
setting  ak  =  b | k j  yields  equation  (5.3). 

Conversely,  Theorem  5.1  implies  that  if  a  satisfies  equation  (5.3)  for  all  polyno¬ 
mials  of  degree  <2N  —  1  and  has  symmetric  support,  then  it  must  be  symmetric. 
Clearly  (5.17)  must  be  also  be  satisfied.  From  (5.2),  this  is  equivalent  to  (5.12a)  and 
(5.14)  with  V2  f(z)  replacing  B(z)  where  f(z)  is  of  degree  2N  —  1  in  z  and  the  also  in 
z  .  Letting  n  =  1  (or  0)  in  (5.14a),  we  see  that  f(z)  has  2N  roots  at  z  =  —1.  Since  f 
is  symmetric,  f(z)  must  also  have  the  form  (5.12b).  We  conclude  that  V2  f  =  h  *  h+ 
where  h  is  a  Daubechies  filter.  Thus, 

Theorem  5.2 

There  is  a  one-to-one  correspondence  between  the  squares  of  the  Daubechies 
orthonormal  waveiet  filters  h  of  length  2N  and  the  Lagrange  ' a  trous  filters  f  of  length 
4N  -  1  given  by  f  =  h  *  ht/V/2. 

Note  that  one  can  compute  the  h  of  length  2N  by  taking  all  possible  square  roots 
of  the  Lagrange  'a  trous  filters  f.10  f  is  easily  computed  from  (5.6)  where  its  even 

^  During  the  revision  of  this  report,  it  was  brought  to  the  author’s  attention  that  an  implicit 
relationship  between  the  squared  filters  and  Lagrange  interpolation  had  been  independently 
noted  in  private  conversations  between  I.  Daubechies  and  Ph.  Tchamitchian 
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components  are  given  by  f2k  =  S^.q/'S/I,  its  odd  positive  components  by  f2ic_i  =■  bk  for  k 
=  0  to  N,  and  for  odd  negative  k  by  symmetry.  Also,  the  spectrum  of  f, 
f(elw)  =  |  h(eIW)  |2,  presents  a  convenient  method  of  computing  the  power  spectra  of 
the  h’s.  In  another  vein,  since  the  h  are  maximally  flat  filters  (i.e.,  have  the  same 
number  of  vanishing  derivatives  at  z  =  1  and  z  =  —1),  Theorem  5.2  shows  that  a  maxi¬ 
mally  flat  filter  is  a  Lagrangian  interpolator,  a  fact  which  may  aid  in  the  design  of 
such  filters.  [14] 
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6.  WAVELET  FILTERS  IN  SIGNAL  PROCESSING 


This  section  has  its  roots  in  a  question  which  originally  motivated  the  author  to 
undertake  this  study:  Inasmuch  as  the  'a  trous  and  Mallat  algorithms  share  the  same 
recursions,  why  not  choose  the  Daubechies  filters  since  they  enjoy  the  additional 
advantage  of  orthonormality?  The  strongest  arguments  in  favor  of  orthonormality 
seem  to  be  mathematical  elegance,  ease  of  inversion,  and,  more  subtly,  good  numeri¬ 
cal  properties.  The  major  drawback  is  a  lack  of  flexibility  in  filter  design,  in  particu¬ 
lar,  an  essentially  fixed  relative  bandwidth.  On  the  other  hand,  we  have  seen  that  the 
DWT  has  a  firm  analytical  basis  independent  of  the  'a  trous  approximation,  even  in 
the  nonorthogonal  case.  In  the  present  section,  we  shall  briefly  examine  the  issues  of 
inversion,  boundedness,  and  adjustable  relative  bandwidth.  (Another  fundamental 
issue,  regularity  of  the  associated  wavelet  functions,  appears  in  appendix  B.)  In  par¬ 
ticular,  we  wish  to  establish  conditions  under  which  these  properties  obtain, 
expressed  directly  in  terms  of  the  filters  f  and  g  rather  than  in  terms,  for  example,  of 
the  scale  function  </>(t)  or  wavelet  t£(t).  This  view  is  in  the  spirit  of  the  DWT  as  an 
entity  in  its  own  right,  and  it  is  certainly  a  necessary  element  in  deciding  which  filters 
to  use  in  practice. 

Any  software  realization  of  the  wavelet  transform  only  implements  a  finite 
number  of  octaves.  Mathematically,  this  reduces  inversion  to  an  algebraic  question, 
one  of  finding  filters  which  satisfy  certain  (not  overly  restrictive)  equations.  How¬ 
ever,  other  considerations  begin  to  come  into  play.  Exact  inversion  requires  finite 
filters,  and,  even  then,  exceedingly  long  filters  may  not  be  useful.  Moreover,  the  con¬ 
strained  problem  is  considerably  more  difficult  to  solve.  An  alternative  approach, 
approximation  by  truncated  infinite  filters,  might  be  acceptable,  but,  once  again,  prac¬ 
tical  considerations  dictate  that  the  filters  decay  quickly.  Similarly,  the  behavior  of 
the  DWT  at  infinity  (i.e.,  w‘  as  i  goes  to  infinity)  becomes  relevant.  For  example,  the 

condition  X  gn  =0,  the  discrete  counterpart  of  (1.2),  is  not  necessary  for  inversion 
n 

of  a  finite  number  of  stages.  However,  it  is  necessary  for  finite  energy  and  bounded¬ 
ness,  which  are  desirable  properties  inasmuch  as  they  reflect  directly  on  the  numeri¬ 
cal  stability  of  the  algorithm  and/or  its  inverse  (see,  for  example,  [3]). 

INVERSION 

To  invert  either  the  decimated  or  undecimated  discrete  wavelet  transform11  it 
suffices  to  invert  a  single  stage  (octave);  that  is,  to  find  s',  given  s1+1  and  w1+1  or  w1+1. 
The  equations  for  inverting  the  decimated  algorithm  are  exactly  analogous  to  those  for 
the  Mallat  algorithm  pictured  in  figure  2.3.  One  seeks  two  filters  p  and  q  which  invert 
a  single  stage  of  the  decimated  DWT  in  figure  1.1;  i.e.,  such  that 

s1  =  PDsi+1  4-  QDsi+1 

=  (PD)  ( AF)  s'  +  (QD)UG)s>  .  (6.1) 


TT 


More  properly  to  invert  a  finite  number  of  stages,  see  next  subsection. 
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Equivalently, 


(PD)(AF)  +  (QD)(AG)  =  I  ,  (6.2) 

where  I  is  the  identity  matrix.  This  type  of  equation,  which  in  the  frequency  domain 
may  be  separated  into  two  equations  comparable  to  (2.14a)  and  (2.14b),  has  been 
treated  extensively  (but  not  exhaustively)  in  the  subband  coding  literature  (cf.,  [12], 
[13],  or  even  [8]).  The  QMF  filters  of  the  Mallat  algorithm  satisfy  (6.2)  with 
g z(w)  =  fz(aH-7r),  p  =ff,  and  q  =  g+  (i.e,  equations  (5.10)  and  (5.11)).  A  less  res¬ 
tricted  class  is  of  the  form  p  =  f*  and  q  =  g*.  The  general  class  of  filters  satisfying 
(6.2),  so-called  biorthogonal  filters,  are  examined  in  [26]  and  [27].  It  should  be 
emphasized  that  for  perfect  reconstruction  in  applications  all  filters  must  be  of  finite 
length.  This  does  not  imply  that  infinite  filters  implemented  by  their  truncations  are 
not  worthy  of  consideration.  [26] 

For  the  undecimated  algorithm  the  requirements  for  inversion  are  much  less 
stringent.  In  order  to  invert  a  stage  of  the  algorithm  of  figure  3.2,  the  filters  p  and  q 
need  only  satisfy  (cf.,  figure  6.1) 

s1  —  p  *  s1+1  +  q  *  wi+1 

=  (p  *  f  4-  q  *  g)  *  s1  .  (6.3) 

That  is, 

p*f  +  q  *  g  =  S  (6.4) 

where  the  Kronecker  delta,  6£60  m,  is  the  identity  for  convolution.  This  is  a  single 
equation,  and  much  more  tractable  than  (6.2).  If  the  polynomials  formed  by  the  z- 
transforms  f(z)  and  g(z)  are  relatively  prime,  one  may  apply  the  Euclidean  algorithm 
for  the  greatest  common  divisor  (in  this  case,  one)  to  find  p  and  q.  It  has  the  advan¬ 
tage  that  finite  f  and  g  lead  to  finite  p  and  q.  Another  method  is  simply  to  solve  the 
equation  in  frequency  space, 

P2(w)  f2(w)  +  q2(w)  gz(w)  =  1  .  (6.5) 

There  is  almost  too  much  flexibility  in  solving  this  equation,  although  it  becomes 
much  more  restrictive  if  one  demands  that  the  filters  be  finite  or  rapidly  decreasing. 
Once  again,  a  popular  choice  [28]  is  f2  f2  4-  g2  jE  =  1,  which,  for  example,  can  be 
solved  for  g2  by  taking  the  square  root  of  1  —  |f“|  as  long  as  |f2(w)|  <1,  (  or,  vice- 
versa,  it  can  be  solved  for  f2).  The  Daubechies  (QMF)  filters  h/V2  and  g/V2  cer¬ 
tainly  satisfy  this  equation  so  that  inversion  for  the  undecimated  version  of  the  Mallat 
algorithm  is  immediate.  Another  case,  useful  in  signal  processing,  is  to  choose  f  to 
be  'a  trous,  p  =  ff/2,  and  g  any  filter  with  nonvanishing  spectrum  except  possibly 
where  |fz(w)|  equals  V2  (cf.,  next  subsection).  Important  questions  of  numerical 
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stability,  filter  lengths,  etc., certainly  remain  to  be  answered,  but  are  well  beyond  the 
scope  of  the  present  report. 


s' 


f 


ni+l 


S 


i+1 


g 


¥ 


w1+1 


q 


w1+1 


Figure  6.1  Illustration  of  a  single  stage  and  its  inverse  for  the 
undecimated  algorithm  found  in  figure  3.1. 


Finally,  before  departing  from  this  subject,  it  should  be  mentioned  that  inversion 
of  the  undecimated  case  in  the  form  of  figure  3.1  also  follows  from  (6.4).  The  invert¬ 
ing  filters  are  just  D'p  and  D'q.  It  is  a  simple  matter  to  verify  that  (6.4)  implies  that 

(D'p)  *  (D'f)  +  (D*q)  *  (D^)  =  b  .  (6.6) 


(Inserting  zeros  in  b  just  yields  b.) 

FINITE  ENERGY  AND  BOUNDEDNESS 

The  discrete  wavelet  transform  is  a  mapping  of  sequences  sn,  n  =  1,  2,  ...  into 
the  space  of  doubly  indexed  sequences  w„,  i,n  =  1,  2,  ...  Finite  energy  for  the  sig¬ 
nal  is  simply 


£  |sn|2  <  oo  .  (6.7) 


One’s  initial  tendency  is  to  look  for  the  same  condition  on  the  wavelet  transform, 
i.e.,  S  lwnl2  <  oo.  However,  the  actual  situation  is  not  quite  so  transparent. 

i,n 

Referring  back  to  the  continuous  case,  we  have  [16] 

/  |s(t)|2  dt  =  /|W(u,6)|2^-  ,  (6.8) 

a~ 

i.e.,  the  weight  (measure)  da  6b  /a1  results  in  units  of  energy.  For  the  DWT,  a  =  2' 
and  d b  is  either  1  (the  undecimated  case)  or  21  (the  decimated  case).  Discretizing  and 
noting  that  da  ja  —  d(lna),  we  see  that 
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da  d b 


(6.9) 


(ln2)  Ai  A b  _  ln2  (—) 
2‘  1  2;  ' 


since  Ai  =  1  .  Observing  that  octave  i  has  a  bandwidth,  up  to  a  constant  factor,  of 
1/21,  we  assign  (6.9)  the  following  physical  interpretation:  The  discrete  wavelet 
transform  is  in  units  of  power/Hz.  Multiplying  by  the  bandwidth  1/21  gives  power, 
and  additional  multiplication  by  the  time  interval  A b  gives  energy.  Note  that  the 
decimated  version  only  outputs  every  2‘  points  so  that  in  that  case  A b  =  21.  In  sum¬ 
mary, 

(i)  |w„|2,  \w'n\2  are  in  units  of  power/Hz. 

(ii)  Octave  i  has  bandwidth  »  1  /2‘. 

(iii)  |w„|2/2'  and  |wj,|2  are  in  units  of  energy. 

One  should  take  care  to  note  that  for  w„,  the  energy  weight  dadb/a 2  is  a  constant 
independent  of  i  so  that  |w„|2  is  discretized  in  a  fashion  so  as  to  be  both  power/Hz 
and  energy. 

Let  1 1  s  1 1 2  4  E  |  sn| 2  be  the  squared  norm  of  s,  and  define 


=  Jim  sup  ~r  ||s'||: 

l— *oo  Z 


(6.10) 


which  corresponds  to  the  DC  energy  (i.e.,  at  u  =  0).  Finally,  define  the  energy  of  the 
DWT  by 


£  =  E  i  ||^|  I2  +  l|s°°||2  (6.11a) 

i  2 

E  =  E  1 1 w* 1 1 2  +  Ms00!!2  .  (6.11b) 

i 

Energy  conservation  takes  the  form  of  the  following  Parseval’s  relationship  for 
discrete  wavelets: 

Definition  6.1 

A  particular  choice  of  filters  f  and  g  is  said  to  be  energy  conserving  if 

||s||2  =  C(E  ^  1 1  w*  1 1 2  +  1 1  s00 1 1 2  )  •  (6.12) 

i  2 

One  may  also  specify  conservation  for  decimated  transforms,  in  which  case  the  2'  is 
dropped. 
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Except  for  some  clarifying  remarks  at  the  end,  we  restrict  the  discussion  in  the 
remainder  of  this  subsection  to  the  undecimated  DWT.  Following  the  above  defini¬ 
tions,  we  see  that  the  wavelet  transform  will  have  finite  energy  if  and  only  if 


£  -^llw'll2  +  I  i  s°°|  |2  <  oo  . 


(6.13) 


This  is  a  necessary  condition  for  the  mapping  w,  from  /2(Z)  to  /2(Z2;  2_I,  1)  to  be 
bounded.  Of  course,  in  practice,  implementations  never  compute  an  infinite  number 
of  octaves.  Nevertheless,  the  property  (6.13)  of  finite  energy  can  be  quite  important. 
Unbounded  transformations  tend  to  have  poor  numerical  behavior  even  when  trun¬ 
cated.  Similarly,  the  inverse  will  not  be  bounded  unless  the  series  (6.13)  is  bounded 
below.  A  wavelet  representation  which  has  these  properties, 


|2  <  E  —Hw1 

I  —  :  01  1  1 


I2  <  B 


(6.14) 


is  called  a  frame  (cf.,  [3],  although  here  the  ambient  Hilbert  space  is  /2  rather  than 
L2(R).)  We  proceed  to  derive  conditions  on  the  filters  f  and  g  for  (6.14)  to  hold. 

Figure  3.1  implies  that 

n  n  n 

/  |s"V)|J  +  /  |w"'(u,)l2  =  /  ( I f, (2^)1 2  +  lt*(2E>)|2)  IsiMI2  .  (6.15) 

—  n  —  7T  —n 

Suppose  that 


max  -(lfz(^’)|2  +  |gz(^’)|2)  <  1 


(6.16) 


Then,  in  the  time  domain,  (6.15)  and  (6.16)  imply 


ill'll2  +  Ill’ll2  <  ilMl2 


(6.17) 


Adding  |  |w’|  |2/2'  to  both  sides  and  repeating  for  decreasing  octaves  implies  that 


S,+'II2  +  E  V  1 1*1 12  <  INI2  • 

i<J  L 


(6.18) 


Finally,  letting  J  go  to  infinity,  we  get  not  only  (6.13),  but  also  the  right  inequality  of 
(6. 14)  with  B  =  1 . 

However,  the  condition  (6.16)  is  much  too  strong,  [’hat  is,  the  transformation 
— *•  (’  g  for  a  large  enough  constant  C  would  cause  (6.16)  to  be  violated  even 
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though  C  has  no  effect  other  than  to  multiply  the  total  energy  by  a  constant.  In  fact 
the  filters  f  and  g  produce  finite  energy  transforms  if  and  only  if  f  and  Cg  yield  finite 
energy.  Thus,  to  have  finite  energy,  it  is  sufficient  to  find  a  C  >  0  such  that 
max  (  |  fz(ce)  | 2  +  C  |  gz(  ca)|2)  <  2.  Such  a  C  exists  provided  that  |fz(oj)|2  <  2  and 

U> 

(1/2)  |gz(w)|2/(l  —  y  |fz(ie)|2)  is  finite;  i.e.,  is  less  than  some  finite  B  =  1/C.  A 
similar  argument  holds  for  the  lower  bound.  If 

min  4(lfz(<-’)|2  +  IgzMI2)  >  1  >  (6-19) 

then  equation  (6.18)  holds  with  the  inequality  reversed  and  the  left  inequality  of  (6.14) 
holds  with  A  =  1.  Once  again,  we  apply  the  trick  with  the  constant  C  and  find  that, 
for  the  inverse  to  be  bounded,  it  is  sufficient  that  there  exist  A  —  1/C  >  0  such  that 

(1/2)  |gzM|2/(l  -  Y  |fzM|2)  >  A.  In  summary, 

Theorem  6.1 

A  sufficient  condition  for  the  undecimated  DFT  w  and  its  inverse  to  satisfy 
(6.14)  (that  is,  to  be  bounded)  is  that,  for  all  u,  |fz(ca)|2  <  2  and 

j  ig^)i2 

0  <  A  <  — - -  <  B  <  oo  .  (6.20) 

1-  jlfzMI2 

To  satisfy  (6.20),  one  must  have  |fz(w)|  —  2  <*=*>  g z(ce)  =0,  and  the  multiplicities  of 
the  corresponding  roots  must  be  identical.  Note,  also,  that  (6.20)  can  be  used  to  give 
an  estimate  of  B/A,  the  so-called  tightness  of  the  frame. 

Whether  these  conditions  are  also  necessary  remains  an  open  question.  One 
can,  however,  show  from  equation  (3.11)  and  an  examination  of  the  power  wz(u;)  at 
ijj  =  0  that  a  necessary  condition  is  gz(0)  =  0  (equivalently,  £  gn  =  0).  This  is  the 

n 

discrete  analogue  of  the  admissibility  condition  (1.2).  The  author  conjectures  that  in 
the  discrete  case  it  is  not  a  sufficient  condition.  (We  remind  the  reader  that  even  in 
those  cases  for  which  the  DWT  is  exactly  the  sampled  WT,  finite  energy  of  the  con¬ 
tinuous  wavelet  transform  does  not  imply  that  of  the  discrete  transform.)  We  do 
have,  however. 

Theorem  6.2 

A  necessary  and  sufficient  condition  for  energy  conservation  (equation  (6.12)j)is 
that  for  all  w 


■7 (I <■*(“••) 1 2  +  jr  I s ,(-')! 2 )  =  1 


(6.21) 
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To  prove  this,  we  may,  without  loss  of  generality,  set  C  =  1  (i.e.,  redefine  g  by  the 
constant  factor  C).  Sufficiency  follows  as  above,  with  inequalities  replaced  by  equali¬ 
ties.  To  prove  necessity,  we  first  note  that  energy  preservation  for  arbitrary  signals 
implies  the  energy  must  be  conserved  for  each  stage.  (For  example,  if  the  signal  is 
s1,  energy  must  be  preserved,  and  since  it  is  preserved  for  s°,  the  first  stage  must 
preserve  energy.)  From  (6.15)  this  implies  (6.21)  with  C  equal  to  one. 

The  decimated  case  seems  to  present  problems.  For  the  above  proofs  to  carry 
over,  the  even  part  of  the  signal  would  need  to  compensate  for  the  lack  of  the  factor 
1 /21  in  (6.11b),  but  |  |Sgven|  |“  7^  1/2  |  | s' |  |2.  This  problem  presents  yet  another  area 
for  additional  research. 

RESOLUTION  AND  RELATIVE  BANDWIDTH 

Considerable  insight  may  be  gained  by  viewing  the  algorithm  in  the  frequency 
domain.  One  stage  of  the  decimated  DWT,  illustrating  equations  (4.11)  from  this 
point  of  view,  is  pictured  in  figure  6.2.  Since  we  are  dealing  with  the  discrete  wavelet 
transforms,  s(z)  =  s(elw)  is  evaluated  on  the  unit  circle.  For  convenience  only  the 
positive  frequencies  are  pictured.  Briefly  the  algorithm  is 

(a)  Bandpass  the  upper  half  of  the  spectrum  to  yield  w1. 

(b)  Lowpass  to  obtain  the  lower  half  of  the  spectrum  (  [0,  7r /2]). 

(c)  Decimate  to  expand  the  lower  half  to  [0,  7r]. 

(d)  Go  to  (a). 

In  somewhat  more  detail:  We  first  obtain  the  high-frequency  information  by  using  g 
to  filter  the  upper  half  of  the  spectrum  of  s'.  The  filter  output  is  w1.  Then,  in  prepara¬ 
tion  for  the  next  octave,  s'  is  lowpass  filtered  by  f.  This  retains  the,  as  yet  unexam¬ 
ined,  low-frequency  contents  and  also  prevents  the  upper  half  of  the  spectrum  from 
aliasing  (i.e.,  contaminating  the  low-frequency  contents)  in  the  dilation  which  follows. 
Finally,  the  operator  A  spreads1-  the  remaining  energy  to  fill  the  spectrum,  producing 
octave  i  +  1.  The  procedure  then  repeats  itself,  s1+1  is  bandpass  filtered  to  get  the 
spectral  contents  at  frequencies  which  are,  in  absolute  units,  one-half  the  frequencies 
of  the  previous  octave. 

A  potential  problem  is  immediately  apparent.  If  the  bandwidth  of  gz(w)  is  less 
than  7t/2,  a  portion  of  the  signal  energy  will  be  discarded;  it  never  appears  in  w1.  One 
possible  remedy  is  to  make  gz  sufficiently  broad;  however,  that  would  limit  the  resolu¬ 
tion.  Alternatively,  we  may  introduce  so-called  voices.  That  is,  we  can  employ  a 
bank  of  filters  of  the  type  g  (see  figure  6.3)  in  order  to  cover  the  entire  upper  half  of 
the  spectrum. 


12 


A, 


which  decimates,  is  a  contraction. 


Its  Fourier  transform  A  is  thus  a  dilation. 
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viewed  from  the  frequency  domain. 


|$(a>)|2 


Figure  6.3  Plot  of  the  power  spectra  of  the  bandpass  filters 
for  four  voices  where  g (oS)  is  given  by  (6.23). 


We  formalize  some  of  these  concepts  using  the  modulated  Gaussian  of  (1.4)  as 
an  example.  With  the  introduction  of  an  additional  parameter  /?,  g(t)  becomes 

g(t)  4  e'"*  e"^ ,2/2  .  (6.22) 

Its  Fourier  transform  is  given  by 

g(w)  =  ^e-(w-">2/2^  .  (6.23) 

We  define  the  bandwidth  of  g(w)  as  twice  the  interval  between  points  for  which  the 
modulus  of  (6.23)  drops  to  1/e  of  its  peak  value,  i.e., 

BW  4  2\/2l3  .  (6.24) 


The  filter  g  is  the  sampled  version  of  (6.22),  that  is, 


Sn 


a  gii/n  g-/?2  nz/2 


(6.25) 


For  convenience,  we  set  the  sample  rate  equal  to  one.  The  following  three 
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restrictions  on  v  and  ft  are  necessary:  First,  in  order  that  g z(u>)  lie  in  the  upper  half  of 
the  spectrum  (cf.,  figure  6.2),  we  require  that 


y  <  u  .  (6.26a) 

Next,  in  order  that  g(t)  be  admissible  and  analytic  (cf.,  [18]),  we  demand 

$<-%;■  (6.26b) 

Under  (6.26b),  gz(a>)  ~  0  for  u  <  0.  ([18]  recommends  ft  <  v/5  as  being  sufficient.) 

Finally,  in  order  that  the  spectrum  not  be  aliased,  we  set 

v  <  7r  —  V2/3  .  (6.26c) 


These  may  be  summarized  in 


max(2n(3,  n/2)  <  u  <  rr  —  V2/9 


(6.27) 


At  this  point  a  word  of  caution  is  advised.  The  bandwidth  of  the  discrete  filter  g 
(e.g.,  (6.25))  is  2\/2  ft  only  when  the  sample  rate  is  1.  Since  the  i,h  octave  is  the 
result  of  i  decimations  by  2,  sampling  the  original  signal  s°  at  a  rate  At  =  1  results  in  a 
Nyquist  frequency  of  n/2'  for  octave  i  (i.e.,  for  s').  Thus,  the  central  frequency  of  g, 
considered  as  a  filter  on  octave  i  is  t//2‘,  and  its  bandwidth  is  2V2/3/2'.  For  this  rea¬ 
son,  it  is  simpler  and  less  ambiguous  to  speak  in  terms  of  a  relative  bandwidth  which 
is  independent  of  sample  rate.  More  precisely,  we  define 


RBW  I 


_ BW _ 

mean  frequency 


(6.28) 


(Any  appropriate  representative  of  the  "center"  of  the  filter  can  replace  the  mean  fre¬ 
quency  in  (6.28).)  In  the  case  of  g,  we  have 


RBW  =  (2 V2/V2‘)  =  2V2 1  > 

(*V2‘)  * 


(6.29) 


Also,  since  7r/2  <  v  <  7r,  we  have  2V2/5//T  <  RBW  < 


4V2/?/tt,  or  approximately 


P  <  RBW  <  2p  . 


(6.30) 


In  view  of  (6.30),  we  shall  consider  the  parameter  ft  as,  essentially,  the  relative 
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bandwidth  of  g. 

The  number  of  voices  M  that  we  would  expect  to  need  to  cover  the  upper  half  of 
the  spectrum  is 


M  ~  i 

The  filter  for  voice  j,  which  we  denote  jV,  is  determined  by  sampling  the  function 
g(t/oJ),  i.e., 

(jv)n  4  g(n/c»j)  for  j  =  0  ,...,  M  —  1  and  a  4  21/M  .  (6.32) 

The  power  spectra  of  the  jV  are  illustrated  in  figure  6.3.  An  alternative  would  be  to 
define  the  voices  as  frequency  translations  of  the  filter  g.  However,  (6.32)  seems  to 
be  a  more  natural  definition  since  it  maintains  the  affine  structure.  (It  is  equivalent  to 
taking  a  nondyadic  value  for  the  dilation  parameter  a  in  equation  (1.1).)  Note  that  the 
bandwidths  of  the  voices  decrease  with  j,  and  the  spectral  spacing  2V2 ft /oj  differs 
somewhat  from  that  assumed  in  (6.31). 

Control  over  the  relative  resolution  is  important  because,  although  the  absolute 
resolution  (0/2')  improves  at  higher  octaves,  at  the  same  time,  the  analysis  frequency 
decreases.  If  one  wishes  to  improve  resolution  at  a  given  frequency,  one  has  to  better 
the  relative  resolution,  i.e.,  decrease  0.  On  the  other  hand,  the  standard  tradeoffs 
apply  to  choosing  0.  Small  0  increases  the  relative  resolution  but  also  requires  more 
voices  and  a  longer  filter  g;  thus,  more  computation.  (The  duration  of  g(t)  is  on  the 
order  of  2\f2/0  so  that  to  provide  a  reasonable  approximation  to  (6.22),  the  length  of 
g  must  be  proportional  to  1  /0.)  Moreover,  the  increased  length  of  g  implies  a  wor¬ 
sening  of  the  (relative)  time  resolution.  The  time-bandwidth  product  is  bounded 
below  by  the  uncertainty  principle,  and  no  amount  of  computation  will  simultaneously 
produce  arbitrarily  small  time  and  frequency  resolution  in  a  single  wavelet  transform. 

With  respect  to  choosing  the  lowpass  filter  f,  we  note  that  the  spectrum  of  a 
longer  filter  f  will  generally  have  a  a  sharper  cutoff.  This  cutoff  is  relevant  because  it 
prevents  the  energy  in  the  upper  half  of  the  spectrum  from  leaking  (aliasing)  into  the 
lower  half  under  the  decimation,  A.  For  most  applications,  a  Lagrange  'a  trous  filter 
of  length  7  (N  =  2)  is  sufficient.  [5]  One  can,  of  course,  also  use  the  asymmetric 
filters  h.  Their  discrimination  of  temporal  direction  seems  intriguing,  but  remains 
uninvestigated. 

Finally,  how  do  the  above  considerations  relate  to  orthonormal  wavelets?  The 
power  spectra  of  filters  h  and  g  obeying  (2.14)  satisfy  (cf.,  [8],  [20]) 

IM-’)l2  +  lgz(*>)l2  =  2  (6.33a) 

and 

I  MO)  |  =  |gz(±7T)|  =  V2  .  (6.33b) 
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It  follows  that,  for  positive  frequencies  (likewise,  for  negative  frequencies),  g  and  h 
must  each  maintain  a  bandwidth  on  the  order  of  7t/2.13  Thus,  in  exchange  for  ortho¬ 
normality  one  relinquishes  control  over  bandwidth.  The  relative  bandwidth  is,  essen¬ 
tially,  fixed  at  (tr/2)  /  (3tt/4)  =  2/3.  On  the  other  hand,  if  one  wishes  to  give  time  and 
frequency  localization  equal  weight  (e.g.,  fi  =  1  so  that  bandwidth  =  duration  =  2V2), 
a  relative  bandwidth  in  the  neighborhood  of  2/3  is  in  a  sense  optimal. 


13  - ,  . 

The  larger  the  value  of  N,  the  more  rapid  the  asymptotic  convergence  of  V'O*')  to  zero;  i.e., 

as  ui  — ►  oo.  [8]  This  hints  at  a  smaller  bandwidth  for  i’  and,  hence,  also  for  g,  for  large  N. 
However,  the  speed  at  which  the  0(tv)  fall  off  near  w  —  0  appears  to  be  fairly  insensitive  to  N. 
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7.  CONCLUSION 


We  have  seen  that  the  'a  trous  algorithm  bears  an  intimate  relationship  to 
Mallat’s  multiresolution  algorithm.  Originally  devised  as  a  computationally  efficient 
implementation,  it  is  more  properly  viewed  as  a  nonorthogonal  multiresolution  algo¬ 
rithm  for  which  the  discrete  wavelet  transform  is  exact.  Moreover,  the  commonly 
used  Lagrange  ' a  trous  filters  are  simply  the  convolutional  squares  of  the  Daubechies 
filters  for  compact  orthonormal  wavelets. 

From  a  broader  viewpoint,  these  two  algorithms  are  instances  of  the  discrete 
wavelet  transform  (DWT),  which,  in  more  conventional  terms,  is  simply  a  filter  bank 
utilizing  decimation  and  two  filters.  There  are  two  basic  versions  of  the  DWT,  one  of 
which  is  simply  the  decimated  output  (octave  i  is  decimated  by  2')  of  the  other.  The 
decimated  DWT  is  characterized  by  octaves:  (a)  obtained  by  alternating  a  lowpass 
filter  f  with  decimation,  and  (b)  tapped  by  a  bandpass  filter  g  to  produce  the  output. 
The  undecimated  DWT  inserts  i  zeros  between  the  elements  of  the  filters  at  octave  i 
in  lieu  of  decimation.  (In  the  case  of  voices,  several  g’s  are  used.)  Finally,  we  note 
that  under  very  general  conditions,  there  exists  a  function  #(t)  such  that  the  fiber 
bank  outputs  w„  correspond  to  :he  sampled  wavelet  transform 
If  —  t 

J  ^(~i "  —  n)  s(0  dt,  thus,  justifying  the  terminology  discrete  wavelet  transform. 

The  personality  of  a  given  DWT  is  distinguished  by  the  choice  of  filters.  If  f 
satisfies  the  'a  trous  condition  f2n  =  <WV2,  then  g  is  the  sampled  version  of  ^(t); 
i.e.,  gn  =  ^(n).  If  finite  length  f  and  g  obey  the  constraints  of  the  multiresolution 
algorithm,  then  the  V2'(0e)(21t  —  n)  are  the  compact  orthonormal  wavelets.  A 
number  of  fundamental  constraints  have  been  discussed.  In  various  combinations 
they  have  a  bearing  on  the  regularity  of  the  wavelet  function,  on  the  energy  in  the 
transform  domain,  and  on  the  boundedness  and  invertibility  of  the  transform.  In  par¬ 
ticular,  we  have  provided  a  set  of  conditions  on  the  filters  sufficient  for  the  transform 
and  its  inverse  to  be  bounded.  The  signal-processing  properties  of  the  discrete 
wavelet  transform  depends  particularly  strongly  on  the  choice  of  g.  The  general  con¬ 
straints  mentioned  above  are  not  very  restrictive  on  g;  however,  there  is  considerably 
less  freedom  in  the  orthonormal  case.  In  particular,  if  orthonormality  is  a  require¬ 
ment,  the  half  bandwidth  of  g  (and,  hence,  the  relative  bandwidth  of  the  wavelet)  is 
no  longer  adjustable.  It  remains  fixed  at  approximately  7r/2. 

Many  topics  remain  for  investigation.  Although  considerable  work  has  been 
done  in  finding  filter  pairs  which  have  a  complementary  set  for  the  inverse  transform 
(cf.,  [8],  [26],  [27]),  it  is  far  from  exhaustive.  The  equivalence  of  maximally  flat  filters 
(with  equal  order  roots  at  0  and  n)  with  Lagrange  'a  trous  filters  as  a  design  tool  is 
perhaps  worthy  of  investigation.  Many  of  our  filter  conditions  on  energy  are  suffi¬ 
cient  but  possibly  not  necessary;  a  tight  set  of  necessary  and  sufficient  conditions  for 
boundedness  would  be  desirable.  Finally,  an  investigation  into  the  quality  of  the 
approximation  of  the  DWT  to  the  sampled  WT  in  the  case  where  it  is  not  exact  could 
be  fruitful.  It  would  perhaps  lend  more  insight  into  the  role  of  the  regularity  of  the 
wavelet  function  in  particular  applications. 
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APPENDIX  A  -  PROOF  OF  LEMMA  3.1  AND  I  FMMA  3.2 


Lemma  3.1 


(AF)U  =  (AF)'0ik_2,B 


(A.l) 


Proof 


For  i  =  1,  we  have 


Then,  by  induction, 


(AF)mk  —  ^2m-k  —  (AF)ok_2m  . 


(A.2) 


(AF)1.,  -  (AFy-'j, 

1 

=  ^  ^2m— j  (AF)I-10  k_2i_ij 

=  e  fH  (AFy-'^-^j-A 

=  S  f-i  (AFy-'i>k.2,m 


(A. 3) 


Lemma  3.2 


E  (AF)U  eil" 


i— 1 

ei2W  n  fz(2»w) 

j-o 


(A-4) 


Proof 


For  i  =  1,  we  have 


E  (AF)nk  eik"  =  E  f2n_k  eik" 

k  k 

=  E  ei2nw  f ’  (w) 

k 


(A. 5) 


Then,  by  induction, 


E  (AFynk  eik 


=  E  (AF)nk  eik(2i_1^  FT  fz(2 )u) 

k  j— 0 

i-1 

=  e,2W  n  fz(2'cj)  . 

j-0 


(A. 6) 


A-l 


APPENDIX  B  -  SUMMARY  OF  FILTER  CONSTRAINTS 


The  discrete  wavelet  transform  w'  and  the  decimated  discrete  wavelet  transform 
w1  (or  d1)  are  defined  for  arbitrary  filters  f  and  g  by 

si+1  =  (D'f)  *  s*  (B.la) 

w1  =  (D'g)  *  s'  ,  (B.lb) 

and 


respectively.  Also, 


si+1  =  A(f  *  s') 
w'  =  g  *  s’ 
d1+1  4  A  w'  , 


(B.2a) 

(B.2b) 

(B.2c) 

(B-3) 


In  essentially  all  applications  f  is  a  lowpass  filter  and  g  is  highpass.  This  rather  vague 
Ratification  is  quantified  below. 

That  is,  in  addition  to  the  above  definitions,  it  is  expedient  to  impose  auxiliary 
conditions  on  the  filters  to  insure  (a)  that  the  DWT  is  related  to  some  WT  with  a  rea¬ 
sonably  behaved  scale  function  <£(t);  (b)  that  the  transform  have  finite  energy  and  be 
a  bounded  transformation;  and,  often,  (c)  that  it  be  invertible.  The  algebraic  condi¬ 
tions  for  invertibility  are  found  in  equations  (6.2)  and  (6.4).  At  the  time  of  this 
report,  no  single  set  of  necessary  and  sufficient  conditions  exist  for  the  satisfaction  of 
(a)  and  (b).  Indeed,  the  definition  of  "reasonable  behavior"  of  the  scale  function  ulti¬ 
mately  depends  on  the  application.  In  an  attempt  to  provide  some  degree  of  organiza¬ 
tion,  we  first  list  the  candidate  constraints,  loosely  labeled  as  lowpass,  highpass,  or 
energy  conditions.  We  then  summarize  their  consequences.  If  either  of  the  filters  is 
infinite  it  is  assumed  to  satisfy  the  decay  condition  (cf.,  [8],  [25]) 

3f  such  that  Y  |fn|nf  <  oo  .  (B.4) 

n 


CANDIDATE  CONSTRAINTS 


(i)  lowpass 


Y  fn  =  V2 

n 

(i.e.,  yfz(0)  =  1) 


(B.5) 


B-l 


(ii)  energy 


y|fz(w)|2  <  l  (B.6) 

(iii)  lowpass 

-feUu)  =  (l+e-)N^(W)  (B.7) 

where  |-7(o;)|  <  C  <  1 
(Note  that  (B.7)  implies  that  fz(7r)  =  0) 

(iv)  energy:  complementary  lowpass/highpass  pair 

J  |gz(^)|2 

0  <  A  <  — =— -  <  B  <  oo  .  (B.8) 

1-jlfzMI2 

Equation  (B.8)  implies 

(a)  highpass 

E  g„  ^  0  (B.9) 

n 

(i.e.,  fe(0)  =  0) 

(b)  lowpass/highpass 

y  I  tz(v) I  =  1  =*>  gz(«)  =  0  (B.lOa) 

gzM  =  0  =>  y|f»|  =  1  (B.lOb) 

(v)  energy 

j(|fz(w)|2  +  c  |gz(w)|2)  =  1  . 


(B.ll) 


IMPLICATIONS 


Necessity  for  pointwise  convergence  of  (4.1)  to  <^(c j):  (B.5) 

Sufficiency  for  (4.1)  and  (4.4)  to  converge  in  LJ(R),  and  L2(R)  to  continuous 
0(t)  and  respectively:  (B.5),  (B.6),  and  (B.7). 

This  is  one  of  the  central  results  of  [8],  which  also  includes  an  examination  of 
the  decay  of  <^>(t)  and  other  regularity  properties.  Note,  that  an  important  class 
of  wavelets  which  do  not  fall  under  the  domain  of  this  theorem  is  the  Haar 
wavelets  (cf., equation  (2.28a)).  Pointwise  convergence  still  holds  for  the  Haar 
wavelets,  but  they  are  not  continuous. 

Necessity  for  for  finite  energy  :  (B.9). 

Sufficiency  for  finite  energy  and  that  the  transformation  be  bounded:  (B.6)  and 
B  <  oo  in  (B.8). 

Sufficiency  for  a  bounded  inverse:  (B.6)  and  A  >  0  in  (B.8). 

Necessity  and  sufficiency  for  energy  conservation  :  (B.ll). 
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